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ABSTRACT 


The  Structural  Synthesis  Transformation  is  used  to  conduct 
structural  system  identification  in  the  frequency  domain.  For 
spatially  complete  cases  where  each  of  the  frequency  response 
functions  at  every  degree  of  freedom  of  each  of  the  coordinates  of 
the  modeled  system  are  available  it  is  shown  that  the  theory  exactly 
identifies  all  modeling  errors.  For  spatially  incomplete  cases  where 
the  frequency  response  functions  are  available  only  at  a  proper 
subset  of  the  degrees  of  freedom  of  the  finite  element  model,  single 
mode  solutions  are  computed  over  intervals  about  the  modes  of  the 
experimental  system  using  matrices  and  complex  valued  line 
integrals.  Methods  of  forming  multiple  mode  solutions  from  the 
single  mode  solutions  are  explored. 
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I.     INTRODUCTION 

The  Finite  Element  (FE)  method  is  a  proven  tool  for  modeling  structural  dynamic 
systems.  As  the  complexity  of  the  system  increases  the  FE  model  may  not  accurately 
reflect  the  dynamic  behavior  of  the  system  To  determine  the  extent  to  which  the  FE 
model  accurately  describes  the  physical  system,  a  comparison  of  the  dynamic  response  or 
modal  parameters  of  the  system  as  predicted  by  the  FE  model  and  the  response  or  modal 
parameters  of  the  physical  system  as  determined  by  measurements  of  the  dynamic 
response  of  the  system  can  be  made.  Such  a  comparison  can  easily  point  out  the 
differences  in  the  dynamic  behavior  of  the  FE  model  and  the  physical  system  but  fail  to 
provide  the  necessary  corrections  to  the  FE  model  that  will  provide  a  more  accurate  FE 
model  prediction  of  the  dynamic  behavior  of  the  physical  system. 

Structural  system  identification  refers  to  procedures  used  to  identify  finite  element 
modeling  errors  using  dynamic  test  data.  Localization  is  the  process  of  identifying  those 
degrees  of  freedom  (DOF)  of  the  FE  model  whose  impedance  differs  from  that  of  the 
physical  system  We  refer  to  this  set  of  DOF  as  the  error  set.  Identification  is  the  process 
of  finding  matrices  AK,  AM,  and  AC  that  are  corrections  to  the  FE  model  stiffness,  mass, 
and  damping  matrices.  When  corrections  are  installed,  the  frequency  response  of  the 
corrected  FE  model  better  predicts  the  frequency  response  of  the  experimental  system, 
and  hence  the  modal  parameters  of  the  corrected  FE  model  better  predicts  those  of  the 
experimental  system. 

Structural  system  identification  can  be  conducted  using  either  modal  or  frequency 
domain  methods.  This  thesis  investigates  a  frequency  domain  method  based  on  the 
structural  synthesis  transformation  (SST)  as  outlined  in  Reference  (1).  When  measured 
dynamic  response  data  is  available  for  each  DOF  of  each  coordinate  of  the  FE  model  the 
identification  is  termed  spatially  complete  and  the  SST  can  be  used  to  exactly  determine 
which  elements  of  the  FE  model  are  in  error.  Additionally,  the  SST  can  be  used  to 
compute  correction  matrices  AM,  AK,  and  AC  that  can  be  used  to  correct  the  mass, 
stiffness,  and  damping  matrices  of  the  FE  model.  The  dynamic  response  of  the  corrected 


FE  model  exactly  matches  that  of  the  physical  system.  In  the  more  frequent  case  that 
dynamic  response  data  is  not  available  at  every  DOF  of  the  FE  model,  the  identification  is 
termed  spatially  incomplete.  In  this  case  the  SST  provides  a  frequency  dependent 
solution. 

SST-based  structural  system  identification  uses  the  frequency  response  function 
(FRF)  of  the  structure  under  consideration.  The  FRF  is  a  quantitative  description  of  the 
dynamic  behavior  of  the  structure.  To  obtain  the  FRF,  known  harmonic  excitation  forces 
are  applied  and  the  resulting  harmonic  response  of  the  structure  measured.  The  ratio  of 
excitation  forces  to  response  at  a  coordinate  evaluated  at  each  frequency  in  a  specified 
bandwidth  defines  the  FRF. 


n.     THEORY 


A.     IMPEDANCE  DESCRIPTION 

The  impedance  model  of  a  given  physical  system  can  be  defined  by  the  relationship 
of  structural  displacement  with  respect  to  an  applied  force, 


IX 


iya         >ya 


(2.1a) 


The  harmonic  force  and  response  vectors  are  denoted  by  "f '  and  "x"  respectively. 
These  vectors  and  the  impedance  matrix,  Z,  are  in  general  complex- valued  and  frequency 
dependent.  Subscripts  "i"  and  "c"  denote  non-error  and  error  DOF  respectively.  The 
superscript  "a"  denotes  quantities  calculated  from  a  FE  (analytic)  model.  If  the  values 
were  obtained  from  experimental  test  data,  the  superscript  would  be  "x".  Thus,  for  the 
experimental  model  the  impedance  relationship  would  be  given  by: 


IX 
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K  cc 


(2.1b) 


In  pratice,  the  elements  of  the  experimental  impedance  matrix  of  Equation  (2.1b) 
are  unmeasurable  quanties.  To  see  this  we  expand  row  j  of  Equation  (2.1b)  to  obtain 


fj=ZjlXl+Z%X2+'--+Z%X„ 


(2.2) 


To  measure  an  element  zxn  of  Z",  requires  that  we  impose  a  unit  displacement  at 

coordinate  xi  while  physically  holding  all  other  other  coordinates  at  zero  displacements. 
This  is  not  physically  possible.  Assuming,  for  purpose  of  definition,  the  availability  of  the 
experimental  impedance  matrix,  the  quantitative  difference  between  the  analytical  and 
experimental  systems,  as  a  function  of  frequency,  is  described  by  the  impedance  error 


matrix.  It  is  defined  by  the  difference  between  the  analytical  and  experimental  impedance 
matrices.  The  error  impedance  matrix  relationship  is  defined  as: 


(2.3) 
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B.     STRUCTURAL  SYNTHESIS  TRANSFORMATION 

Since  the  experimental  impedance  matrix,  T,  is  in  general  unavailable,  frequency 
domain  structural  synthesis  is  used  to  identify  the  impedance  error  matrix  using  FRF  data 
exclusively.  A  structural  synthesis  transformation  is  constructed  from  AZ  of  Equation  (2.3) 
which  encompasses  the  FE  model  errors.  This  transformation  is  applied  to  the  finite 
element  model  to  produce  an  experimental  system  FRF. 

The  FRF  relates  structural  response  to  applied  excitation.  Given  a  FE  model  with 

impedance  matrix,  Z3,  the  FRF,  H°,  is  the  matrix  inverse  of  the  impedance  matrix  2°  of 

equation  (2.1a).  For  a  static  system  (£2=0)  the  FRF  is  simply  the  flexibility  matrix  (inverse 

stiffness  matrix).  Using  the  notation  of  Equation  (2.1)  we  may  partition  H°  as  follows: 
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(2.4) 


Generally,  the  "c"  response  coordinates  experience  applied  forces  due  to  both  error 
impedances  and  externally  applied  forces,  whereas  "i"  response  coordinates  experience 
only  externally  applied  forces,  such  that, 
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AZ 


(2.5a) 


and 


f,=ft 


(2.5b) 


Expanding  Equation  (2.4)  and  substituting  the  relations  of  Equation  (2.5)  yeilds 
the  following  relationships: 


a  rext    ,     rj-a  rext    .     rja  /-AZ 


Xi  ~  niiJi       +  MicJc       +  nicJc 


xc = K/r + kjt + k/c 


(2.6a) 


(2.6b) 


If  we  include  a  copy  of  Equation  (2.6b),  in  expanded  matrix  notation,  Equation 
(2.6)  reflects  the  three  harmonic  excitation  terms  to  be  considered,  i.e., 
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(2.7) 


Response  coordinates  "c"  and  "i",  which  are  due  to  external  forces,  will  hereafter 
be  referred  to  as  "e"  coordinates,  denoting  their  dependence  on  external  force  excitation. 
Consequently,  the  three  excitation  forces  are  condensed  into  two  under  the  identity: 


\f-H/r]  [fr\\ 


(2.8a) 


(2.8b) 


and  Equation  (2.7)  reduces  to 


\Xe 


ce       ncc 


(2.9) 


Equation  (2.3)  shows  that  the  impedance  error  is  defined  as  the  difference  between 
the  analytic  and  experimental  impedance  models.  Hence,  a  transformation  is  required 
which  uses  the  FRF  relationship  of  Equation  (2.9)  to  generate  a  similar  relationship  for  the 


experimental  system.   The  impedance   error  AZ  provides   the   basis   by  which  this 
transformation  is  developed. 


The  impedance  error  matrix,  AZ(Q),  is  a  function  of  frequency  and  satisfies: 


{/c}  =  -[AZ(Q)]{xc} 


(2.10) 


where 


[AZ(Q)]  =  [[AA:]  +  yQ[AC]  -  Q2[AM]] 


(2.11) 


for  Q  the  forcing  frequency,  /  =  V-T  and  AK,  AC,  and  AM  stiffness,  damping  and  mass 
error  matrices  comparable  to  those  of  the  finite  element  formulation.  The  minus  sign  in 
Equation  (2.10)  reflects  that  the  reaction  forces  imposed  by  impedance  errors  on  the 
baseline  model  are  being  considered.  Substituting  the  relationship 
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into  Equation  (2.9)  yields: 
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simplifying  we  get: 
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(2.14) 


Expanding  Equation  (2.14)  into  two  equations  and  using   "*"  to  denote  a 
synthesized  modified  response  results  in 


xc  =  Hacefe  -  H^AZcl 


(2.15a) 
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(2.15b) 


Rearranging  Equation  (2. 15b)  produces 

[l  +  HaccAZ]x'c=Hacefe  (2.16a) 

Using  the  property  of  the  frequency  response  function, 

xc  =  Hacefe  (2.16b) 

[l+HaccAZyc=xc  (2.16c) 

xc  =[l  +  HaccAZ}~1 lxc  (2.16d) 

Introducing  Equation  (2.16)  into  Equation  (2.15b)  results  in 

*c*  =  Kf€  -  H:CAZ[I  +  KAZ}-1  xc  (2.17a) 

xc  =  H°eefe  -  HaecAZ[l+  HaccbZ\X Hacefe  (2.17b) 

Once  again  we  recall  the  property  of  a  FRF, 

*e*  =  Kefe  (2.18a) 

Combining  Equations  (2.17b)  and  (2.18a)  yields 

K  =K  -HaecAZ[l+H:cAZ]~lHace  (2.18b) 

Noting  that 

[/+  i/^AZ]"1  =  [(AZ-1  +  H^AZ]'1  (2.19a) 

and  applying  the  matrix  property 

([<#])- =[*]-'["]-'  (2.19b) 

we  get  that 

lTee  =  Haee-  //^AZ"1  +  K]'1^  (2.20a) 


Replacing  the   superscript   "*",   which  denotes  the   structures's   synthesized  coupled 
response,  with  the  superscript  "x"  to  indicate  the  test  system  response  we  arrive  at 

/£  =  /£-  h:c[aZ~1  +  H°cc\lHace  (2.20b) 


In  full  matrix  notation  we  have 


'HI 

K 

~m 

Hi 

~K 

Hi 

Hzcc_ 

Haa 

Hacc_ 

K 

[AZ-  +  H°X 


~\T 


(2.20c) 


Equation  (2.20b)  is  the  structural  synthesis  transformation  equation.  When  the 
experimental  system  FRF,  If  ,  is  available  the  SST  can  be  used  to  identify  a  frequency 
dependent  impedance  error  matrix  [AZ(Q)].  Additionally,  using  Equation  (2.9),  [AZ(Q)] 
can  be  decomposed  into  constitutent  stiflhess,  mass,  and  damping  errors. 

C.      FREQUENCY  DOMAIN  LOCALIZATION 

We  may  rewrite  Equation  (2.20b)  as 


where 


and 


A#«  =  H°ecD-XH°ce 


&Hee   ~  H°ee  ~  H^ 


Z)  =  [AZ-1+//;e] 
We  define  the  localization  matrix  L  as 


L  -  Zaee-  AHee  •  Z°e 


using  the  expression  of  Equation  (2.22a)  in  Equation  (2.23)  we  can  rewrite  L  as 


L  =  Z:e.HaecD~lHace-Z:e 


Expanding  the  "e"  coordinate  set  into  error  and  non  error  coordinates  we  get 


(2.21) 


(2.22a) 


(2.22b) 


(2.23) 


(2.24a) 
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L  = 


rya  'yd 


E>i 


hi  h: 


ya        ya 


(2.24b) 


Noting  the  the  frequency  response  matrix  is  the  inverse  of  the  impedance  matrix,  i.e., 


za 

za 

~H° 

H°c 

'i 

0" 

za 

za 

_H« 

K_ 

0 

I 

(2.25a) 


all  mixed  product  coordinates  of  Equation  (2.24b)  must  be  zero  and  L  simplifies  to 

"0      0 
0    D~x 


L  = 


@Q=Qi 


(2.25b) 


D.     ERROR  IMPEDANCE 

We  can  solve  Equation  (2.20b)  for  the  impedance  error  [AZJ.  From  Equation 
(2.25b)  terms  of  Equation  (2.20b)  associated  with  non  error  coordinates  may  be  assumed 
to  be  zero.  We  get  the  following  form  of  the  impedance  error  matrix 


-i 


where 


[az]  =  ([K]- [*£])'   @Q=Qi 


^c]=([//:]-1[A^|/f:]-1)"1 


(2.26a) 


(2.26b) 


Let  S={Qi,  Q2,  —  ,  0.D}  be  a  set  of  frequencies  where  Qi  <  Q2<...<^n-i<^n.  If 
for  each  i=  2,3, ...,  n-1,  we  apply  Equation  (2.10)  at  each  of  the  frequencies  Qi_i,  Qj ,  and 
Qi+i  and  assemble  the  resulting  three  equations  into  a  system  of  three  equations  in  three 
unknowns,  we  get  the  matrix  equation 


AZC(Q,) 
AZC(Q,+1) 


I    -«f+17    jQmI 


AKC 
AMC 
AC 


(2.27) 


Equation  (2.27)  can  be  used  to  decompose  the  frequency  dependent  impedance 
error  into  constitutent  strShess,  mass,  and  damping  error  matrices,  AK,  AM,  and  AC  at  the 
frequencies  Qi,  i=2,3,...,n-l.  These  constituents  matrices  are  in  general  frequency 


dependent.  Denoting  the  solution  of  Equation  (2.27)  by 


A£C(Q,) 
AMC(Q,) 
ACC(Q,) 


for  i=2,3,...,n-l  we 


obtain  for  each  Cl[  i=2,3,...,n-l,  error  stifihess,  mass,  and  damping  matrices  AKc(Qi), 
AMc(Qi),  and  ACc(Oi)  such  that 


AZC(Q,)  =  A*C(Q,)  -  Q?AMC(Q,)  +  /Q,AQ(Q,) 


(2.28) 


The  matrices  AKc(Di),  AMc(Oi),  and  ACc(Oi)  can  be  used  to  numerically  correct 
the  stiffness,  mass,  and  damping  matrices  of  the  FE  model  at  the  frequencies  Q; 
i=2,3,...,n-l  in  that  the  FRF  of  the  corrected  FE  model  at  the  frequencies  Qj  approximates 
the  experimental  system  FRF  at  Qj.  To  express  this  symbolically,  if -A",  Af,  and  C  are  the 
stifihess,  mass,  and  damping  matrices  of  the  FE  model  and  If  is  the  FRF  matrix  of  the 
experimental  system  at  Q, 

Hx(n^(Ka  +  AK^-Q^M"  +  AMc(Q;))  +  /Q,(Cfl  +  ACC(Q,)))_1      (2.29) 

where  the  "c"  subscripted  matrices  are  added  at  the  corresponding  error  set  coordinates  of 
the  "a"  superscripted  matrices. 

For  a  general  set  of  frequencies,  E={Qi,  Q2, ... ,  Qn},  we  can  form  the  system  of  n 
equations  in  three  unknowns  given  by: 

AZc(Qj 

AZC(Q,)         /    -Q]l    jX1J\AMc 

AZC(Q„) 


/ 

-a]i 

jQ/ 

~AKC' 

/ 

-Q]I 

JQJ 

AMC 

/ 

-&J 

7*V 

The  solution 


~AKC~ 

-AKc(~y 

AMC 

= 

AMC(3) 

-ACc_ 

ACC(H) 

of  Equation  (2.30)  represents  error  stifihess, 


mass,  and  damping  matrices  which  best  corrects  the  FE  model  in  a  least  squares  sense. 
Equations  (2.27)  and  (2.30)  are  fundamental  to  all  that  follows. 
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m.     SPATIALLY  COMPLETE  STRUCTURAL  IDENTIFICATION 

To  illustrate  the  principles  of  SST  based  frequency  domain  structural  identification 
we  will  use  the  FE  model  of  a  free-free  beam.  To  simulate  the  experimental  system  we 
impose  a  25%  addition  to  the  mass  and  stiffness  of  elements  3  and  4  of  a  10  element  FE 
model.  Figure  3-1  shows  the  finite  element  model  of  the  beam  and  the  spatially  complete 
experimental  system  that  results  from  imposing  the  mass  and  stiflhess  additions  at  element 
3  and  4  of  the  FE  model.  Table  3-1  shows  the  system  frequencies  of  the  analytic  and 
experimental  systems. 
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EXPERIMENTAL  SYSTEM 
MASS  OF  ELEMENT  3  =  1.954X10^6  LBM 
MASS  OF  ELEMENT  4  =1.954X10»-6  LBM 
STIFFNESS  OF  ELEMENT  3  =  417.88  LBF/IN 
STIFFNESS  OF  ELEMENT  4  =  417.88  LBF/IN 
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o    a    o 


a    a    o 
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ANALYTIC  SYSTEM 
NUMBER  OF  ELEMENTS  =  10 
LENGTH  =  60.625 
El  =  74500  LBF-INA2 
SECTIONAL  AREA  =  0.8359  IN*2 
WEIGHT  DENSITY  =  0.05  LBF/IN»3 
ELEMENTAL  MASS  =  1.563X10*-6  LBM 
ELEMENTAL  STIFFNESS  =  334.3  LBF/IN 


f      ^    ROTATIONAL  DOF  AT  NODE 

TRANSLATIONAL  DOF  AT  NODE 


Figure  3- 1  Spatially  complete  Analytic  and  Experimental  Systems. 
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MODE 

ANALYTIC  (Hz) 

EXPERIMENTAL 
(Hz) 

1 

25.41 

24.42 

2 

70.06 

68.41 

3 

137.4 

130.2 

4 

227.5 

216.3 

5 

340.8 

329.4 

6 

478.0 

456.5 

7 

639.9 

612.0 

8 

826.5 

802.6 

9 

1026 

969 

10 

1363 

1323 

11 

1641 

1573 

12 

1981 

1920 

13 

2381 

2283 

14 

2850 

2754 

15 

3397 

3269 

16 

4028 

3912 

17 

4720 

4605 

18 

5377 

5156 

19 

6795 

6791 

20 

6808 

6801 

Table  3-1  System  Frequencies  of  spatially  complete  Analytic  and  Experimental  systems. 

We  will  denote  by  Ma ,  Ka ,  and  Ca  the  mass  stiffness  and  damping  matrices  of 
the  FE  model  and  byM*,  Kx ,  and  Cx  the  mass  stiffness  and  damping  matrices  of  the 
experimental  system.  The  impedance  matrices  of  the  analytic  and  simulated  experimental 
systems  are  given  by: 
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Za(Q)  =  Ka  +jQCa-Q2Ma  (3.1a) 

ZX(Q)  =  Kx  +jQCx  -Q2M*  (3.1b) 

The  FRF  matrices  of  the  analytic  and  simulated  experimental  systems  are  given  by: 

Ha(Q)  =  Za(Q)~l  (3.2a) 

HX(Q)  =  Zx(Q)~l  (3.2b) 

Figure  3-2  shows  a  comparison  of  the  driving  point  FRF  at  DOF  1  of  the  analytic 
and  experimental  system. 
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Figure  3-  2  Analytic  FRF  vs  simulated  Experimental  FRF. 
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-  For  a  given  frequency  Q0,  using  Equation  (2.24a),  we  can  form  the  localization 
matrix  L.  Plotting  the  diagonal  elements  of  L  versus  the  associated  DOF  we  obtain  the 
plot  shown  in  Figure  3.3.  For  each  frequency,  Q,  in  a  given  frequency  range,  we  can 
compute  the  diagonal  of  L  at  Q  .  Assembling  the  diagonals  over  the  frequency  range  of 
our  system  into  a  rectangular  matrix  and  performing  a  MATLAB  mesh  plot  of  the 
resulting  rectangular  matrix  we  obtain  the  surface  plot  shown  in  figure  3.4.  Figure  3.4 
shows  the  frequency  dependency  of  the  localization  matrix  diagonals.  As  our  system  is 
spatially  complete  Equation  (2.25)  forces  all  non  error  coordinates  to  be  zero.  From 
Figure  3-3  we  can  determine  that  the  locations  of  the  nonzero  diagonal  values  are  DOF  9, 
10,  11,  and  12.  We  denote  by  C^  the  set  of  DOF  for  which  L(i,i)  is  non  zero.  For  our 

beam  system  C^  =  {9,10,11,12}  .  Figure  3.5  shows  the  frequency  dependence  of  typical 
error  and  non  error  set  diagonal  elements  of  the  localization  matrix  L. 

Using  the  set,  Cerr,  which  results  from  the  localization,  we  can  perform  a 
partitioning  of  the  FRF  matrix  as  described  by  Equation(2.4).  We  can  now  apply 
Equation  (2.26)  to  compute  AZ  as  a  matrix  function  of  frequency  over  the  frequency 
range  of  our  system.  We  then  use  Equation  (2.27)  to  decompose  AZ  into  its  constituent 
components  AK,  AM,  and  AC.  We  get  exact  solutions  of  error  stifihess,  mass,  and 
damping  as  shown  in  figures  3-7,  3-8  and  3-9.  The  MATLAB  Routine  SST.M  of 
Appendix  A  can  be  used  to  accomplish  the  above  steps. 

For  each  Q  in  the  frequency  range  of  our  system  we  can  form  the  sum: 

Zcorr(Q)  =  Zacc(Q)  +  AZ(Q)  (3.3) 

We  refer  to 

Hcorr={ZcorryX  (3.4) 


as  the  corrected  FRF  of  the  analytic  model.  i/corr  is  the  FRF  of  the  model  that  results  from 
installing  the  corrections  as  identified  by  the  Equation  (2.26a).  Figure  3-10  shows  a 
comparison  plot  of  the  FRF  of  the  corrected  model  and  the  FRF  of  the  experimental 
system.  Figure  3-10  clearly  shows  the  exactness  of  the  SST  solution  in  the  case  of  a 
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spatially  complete  system  The  experimental  and  corrected  model  FRFs  are  identical  to 
within  plot  resolution. 

At  each  frequency,  Q,  in  the  frequency  range  of  our  system  we  can  form  the  sum 

ZZ?(C1)  =  Z°cc{0)  +  AK(Q)  +  jQAC(Q)  -  Q2  AM(Q)  (3.5) 

We  refer  to: 

Hcr={ZccZ'Y  (3.6) 


as  the  constituent  corrected  FRF.  Figure  3-11  is  a  comparison  plot  of  the  constitutent 
corrected  FRF  and  experimental  for  our  system.  The  offset  between  the  two  plots  is 
exactly  equal  to  the  sampling  frequency  AQ. 
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Figure  3-  3  Spatially  complete  localization  matrix  diagonal  at  £1=352  Hz. 
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Figure  3-  4  Frequency  dependence  of  spatially  complete  localization  matrix  diagonals. 
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Figure  3-  5    Frequency  dependence  of  spatially  complete  localization  matrix  error  set 
DOF.  Units  are  Min  (Log  10  of). 
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Figure  3-  6  Spatially  complete  Analytic  Impedance  vs  Experimental  impedance. 
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Figure  3-  7   Computed  stiffness  vs  true  stiffness  for  spatially  complete  beam.  Plots  are 
identical  within  plot  resolution. 
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Figure  3-  8  Computed  mass  vs  true  mass  for  spatially  complete  beam.  Plots  are  identical 
within  plot  resolution. 
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Figure  3-  9   Computed  damping  vs  true  damping  for  spatially  complete  beam.  Plots  are 
identical  within  plot  resolution. 
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Figure  3-  10    Spatially  complete  Experimental  FRF  vs  AZ  corrected  FRF.  Plots  are 
identical  within  plot  resolution. 
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Figure  3-11  Spatially  complete  experimental  FRF  vs  AK,  AM,  and  AC  corrected  FRF. 
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IV.     SPATIALLY  INCOMPLETE  STRUCTURAL 
IDENTIFICATION 

A.     GENERAL  DESCRIPTION 

To  illustrate  SST  based  frequency  domain  structural  identification  when  the 
physical  system  under  consideration  is  spatially  incomplete  we  will  again  use  the  FE  model 
of  a  free-free  beam  We  simulate  the  experimental  system  by  imposing  a  25%  addition  to 
the  mass  and  stifihess  of  elements  3  and  4  of  a  10  element  FE  model.  Figure  4-1  shows  the 
FE  modeled  beam  and  the  spatially  incomplete  system  that  results  if  FRF  data  is  available 
only  at  the  displacement  DOF  of  the  FE  system. 
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ANALYTIC  SYSTEM 
NUMBER  OF  ELEMENTS  =  10 
LENGTH  =  60.625 
El  =  74500  LBF-INA2 
SECTIONAL  AREA  =  0.8359  INA2 
WEIGHT  DENSITY  =  0.05  LBF/IN*3 
ELEMENTAL  MASS  =  1 .563X1 0*-6  LBM 
ELEMENTAL  STIFFNESS  ■  334.3  LBF/IN 
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Figure  4- 1  Analytic  and  experimental  spatially  incomplete  systems. 
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To  obtain  a  simulated  FRF  for  our  spatially  incomplete  beam  we  denote  by 
Ma ,  Ka ,  and  Ca  the  mass  stiflhess  and  damping  matrices  of  the  FE  model  and  by  Mx , 
K* ,  and  Cx  the  mass  stiffness  and  damping  matrices  of  the  experimental  system  Mx , 
Kx ,  and  Cx  are  obtained  by  imposing  15%  mass  and  stiffness  and  a  15%  mass  additions 
to  the  elemental  matrices  of  element  5  and  6  respectively  of  the  FE  model.  The  impedance 
matrix  of  the  simulated  spatially  complete  beam  is  given  by 


ZX(Q)  =  Kx  +jQCx  -  Q2MX  (4-1) 

The  FRF  matrix  of  the  simulated  spatially  complete  system  is  given  by 

HX(Q)  =  Zx(Q)~l .  (4-2) 


We  introduce  the  terminology  Analysis  set  and  Omitted  set  where  the  Analysis  set 
(A-set)  is  that  set  of  DOF  for  which  experimental  FRF  data  is  available  and  the  Omitted 
set  (O-set)  is  that  set  of  DOF  of  the  experimental  system  for  which  experimental  FRF  data 
is  unavailable.  For  our  simulated  experimental  system  the  A-set  consist  of  the  odd 
numbered  translational  DOF  and  the  O-set  consists  of  the  even  numbered  rotational  DOF, 
i.e., 

A -set  =  {1,3,5,..., 21}  (4.3a) 

O-set  =  {2,4,6,..., 22}  (4.3b) 

We  obtain  the  simulated  FRF  of  our  spatially  incomplete  beam  by  physically 
extracting  the  rows  and  columns  of  the  simulated  spatially  complete  matrix,  If,  for  which 
FRF  data  would  be  available.  In  our  system  these  are  the  rotational  DOF  and  all  even 
numbered  rows  and  columns  are  omitted  from  the  simulated  spatially  complete  FRF  to 

obtain  a  spatially  incomplete  FRF  which  we  will  denote  by  Hx  . 

For  a  fixed  Q0,  HX(Q.0)  is  a  square  matrix  of  size  (length^ -set))  by  (length^  - 

set)).  For  our  FE  model  as  currently  defined,  If,  is  of  size  (number  of  DOF)  by  (number 

of  DOF)  and  is,  as  is  true  of  most  real  world  cases,  of  larger  size  than  H X  .  In  order  to 
employ  the  structural  synthesis  transformation  we  need  to  reduce  the  size  of  If  to  that  of 
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Hx  .  To  this  end  we  will  consider  two  reduction  methods,  FRF  matrix  extraction  [Ref.  1], 
and  the  Improved  Reduced  System  as  given  in  [Ref.  2]. 

B.     EXTRACTION  REDUCTION  METHOD 

In  order  to  reduce  H°  by  the  extraction  method  we  simply  extract  from  the  frill 
order  H  those  rows  and  column  which  correspond  to  A-set  coordinates.  Partitioning  the 
impedance  and  full  FRF  matrices  of  our  analytical  system 


Z  = 


H  = 


za 

Za 

za 

_     oa 

Za 

oo_ 

HI, 

Hao 

_Ho\ 

H$o_ 

(4.4a) 


(4.4b) 


and  using  the  identity 


ZH  = 


Za 


#« 

Haao 

'I 

0" 

K 

H°oo_ 

0 

/ 

(4.5) 


we  obtain  the  relationship 


•"oa   ~\^aa        ^ao^oo^oa) 


(4.6) 


We  denote  the  reduced  Analytic  FRF  of  Equation  (4.6)  by  H"  . 

C.      O-SET  SYSTEM 

Equation  (4.6)  relates  the  extracted  reduced  order  FRF  of  the  analytic  model  to  the 
impedance  of  the  full  order  order  model.  Taking  advantage  of  the  identity 


[Z^Y  =  (det[z4~  adj[Z00\  (4.7) 

and   replacing   Z00   by  Kco+jQCoo-jtfMoo   (where   det(«)    and   adj(»)   represent   the 

determinant  and  adjount  respectively),  we  see  that  an  element  H°(Q)  is  large  for  those 

frequencies  Q0  where  Q0  is  an  eigenvalue  of  the  O-set  system,  the  O-set  system  being  that 
FE  model  having  stiffness,  mass  and  damping  matrices  Koo,  Moo,  and  Coo  respectively. 
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D.     IMPROVED  REDUCTION  SYSTEM 


To  use  the  improved  reduction  method  (IRS)  we  first  partition  T  using  the  A  and 
O  sets,  then  adjust  Equation  (2.1)  to  reflect  this  new  coordinate  system  obtaining 


I/.. 


ZQ  rjQ 


(4.8) 


Expanding  Equation  (4.8)  into  two  equations  yields 

f  =  Za  x  +Za  x 


f    =  Z"  X  +  Z"  X 


(4.9a) 


(4.9b) 


O-set  coordinates  are  not  associated  with  FRF  data  measurement  locations  on  the 
physical  structure,  therefore  the  forcing  function  at  O-set  coordinates  can  be  set  to  zero. 
Making  these  substitution  into  Equation  (4.9b)  and  solving  for  the  generalized  structural 
response  coordinates  leads  to: 


x   =  -Z~xZa  x 


I*- 


/ 


M 


(4.10a) 


(4.10b) 


Substituting  these  results  into  Equation  (4.8)  yields, 


10 


yaa        <yo 
ya  rya 


yd       *  yd 

~^oo     ^oa 


{*.} 


(4.11) 


hence 


\Ja  ]  ~  [  ^aa        ^ao^oo  ^oa  j  \Xa  ) 


(4.12) 


When  Q=0  the  Equation  (4.10a)  yields  the  static  reduction  relationship  between 
omitted  and  retained  coordinates  and  is  given  by 


to  =  [-*;X,]to 


(4.13) 
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The  IRS  relationship  is  given  by 

{*.}  =  [~^X  +  TM^Km){xe)  (4.14) 

where 

T  =  K^  -  K'JM^K-JK^  (4.15) 

and  Kstat  and  M^t  are  the  statically  reduced  [Ref.  2,  3]  stiflhess  and  mass  matrices. 

Unlike  the  spatially  complete  case  where  we  only  had  to  consider  two  system  (the 
analytic  and  the  experimental  systems),  in  the  case  of  a  spatially  incomplete  system  there 
are  actually  five  systems  with  which  we  must  concern  ourselves  [Ref  4]:  the  analytic 
system,  the  experimental  system,  the  reduced  analytic  system  that  results  from  conducting 
dynamic  reduction  on  the  mass  and  stiflhess  matrices  of  the  analytic  system,  and  the 
omitted  systems  of  both  the  analytic  and  the  experimental  system.  Table  4-1  shows  the 
frequencies  of  each  the  first  four  of  these  systems. 

Figure  4-2  shows  a  comparision  of  the  analytic  and  experimental  FRF  of  our 
spatially  incomplete  beam.  We  see  that  by  using  IRS  reduction  of  the  analytic  system  those 
modes  above  approximately  1000  Hz  i.e.,  those  modes  associated  with  the  reduced  out 
rotations,  are  not  present.  Figure  4-3  shows  a  similar  comparision  where  we  have  used 
extraction  reduction  and  the  higher  modes  are  present.  In  all  that  follows  we  shall  use  IRS 
reduction  of  our  analytical  systems.  Figure  4-4  shows  the  localization  matrix  diagonal  at 
Q=  196.1  Hz  while  Figure  4-5  shows  the  frequency  dependency  of  the  localization 
diagonals  over  the  frequency  range  our  our  spatially  incomplete  beam.  For  a  spatially 
incomplete  system  the  determination  of  the  locations  of  the  error  coordinates  is  not  a 
clearly  defined  task.  We  shall  not  discuss  the  problem  of  actually  detennining  the  exact 
error  set  in  the  spatially  incomplete  case  and  will  use  our  knowledge  of  the  true  location  of 
the  error  coordinates  to  aid  in  our  localization.  As  we  shall  make  use  of  the  concept  of  the 
size  of  the  error  set  again,  we  will  simply  note  that  using  a  reduction  method  like  IRS 
causes  errors  to  be  'smeared'  in  the  reduced  analytic  model.  Table  4-2  shows  the  A-set  and 
O-set  of  our  simulated  beam  along  with  the  locations  of  the  true  errors  as  well  as  the 
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computed  C-set  at  Q=196  Hz.  The  computed  C-set  is  the  set  of  all  DOF  having  a  diagonal 
entry  whose  absolute  value  exceeds  a  given  tolerance. 


Mode 

Anal  (Hz) 

Exp  (Hz) 

Reduced 
(Hz) 

O  set  (Hz) 

1 

25.41 

25.83 

25.41 

1244.28 

2 

70.07 

70.85 

70.07 

1287.94 

3 

137.45 

137.93 

137.45 

1417.14 

4 

227.55 

227.37 

227.55 

1629.59 

5 

340.84 

341.14 

340.84 

1928.6 

6 

478.0 

477.97 

478.0 

2327.83 

7 

639.9 

639.92 

6403 

2853.66 

8 

826.5 

826.24 

829.7 

3541 

9 

1026.8 

1025.99 

1029.74 

4399.32 

10 

1363.62 

1364.51 

5280.86 

11 

1640.62 

1640.71 

5701.99 

12 

1981.19 

1984.61 

13 

2381.4 

2380.32 

14 

2850.01 

2847.38 

15 

3397.39 

3397.74 

26 

4027.82 

4026.28 

17 

4720.15 

4738.64 

18 

5376.83 

5380.09 

19 

6794.91 

6787.09 

20 

6807.52 

6804.12 

Table  4-1  Analytic,  Experimental,  Reduced,  and  Omitted  System  Frequencies  of  Spatially 
Incomplete  Beam. 
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A  Set 

1 

3 

5 

7 

9 

11 

13 

15 

17 

19 

21 

OSet 

2 

4 

6 

8 

10 

12 

14 

16 

18 

20 

22 

True  C  Set 

5 

6 

7 

8 

9 

10 

Computed      C 
Set 

1 

3 

5 

7 

9 

11 

13 

Table  4-2  Analytic  set,  Omitted  set,  Computed  C  set,  and  True  C  set  DOF  for  spatially 
incomplete  beam. 
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Figure  4-  2  Analytic  FRF  vs  experimental  FRF  for  spatially  incomplete  beam  using  IRS 
reduction. 
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Figure  4-  3    Analytic  FRF  vs  experimental  FRF  for  spatially  incomplete  beam  using 
extraction  reduction. 
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Figure  4-  4  Localization  matrix  diagonal  at  Q= 196.1  Hz  for  spatially  incomplete  beam. 
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Figure  4-  5      Frequency  dependence  of  localization  matrix  diagonals  for  spatially 
incomplete  beam. 
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Figure  4-  6     Frequency  dependence  of  error  and  non  error  set  localization  matrix 
diagonals  for  spatially  incomplete  beam.  Units  are  lbffin  (Log  10  of). 
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Figure  4-  7  Analytic  vs  experimental  impedance  for  spatially  incomplete  beam. 
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Figure  4-  8  Computed  Stiffness  vs  True  Stiffness  for  spatially  incomplete  beam. 
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Figure  4-  9  Computed  Mass  vs  True  Mass  for  spatially  incomplete  beam. 
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Figure  4- 10  Computed  Damping  vs  True  Damping  for  spatially  incomplete  beam. 


41 


0 


-1  - 


-2 


e-3 

o 

£-4 


6. 


-6- 


-7- 


i 

r~ 

i 

I7Y1H7I 

OMENTAL  F  IF 
ECTEDFRF 

CORR 

1 

1 

i 

! 

i 

v 

l\/\ 

i 

i 

i 

200 


400 


600 
Omega  (Hz) 


800 


1000 


1200 


Figure  4- 11  AZ  Corrected  FRF  vs  experimental  FRF  for  spatially  incomplete  beam.  Plots 
are  identical  to  within  plot  resolution. 
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Figure  4-  12    AK/AM/AC  Corrected  FRF  vs  experimental  FRF  for  spatially  incomplete 
beam.  The  offset  of  the  two  plots  is  equal  to  the  sampling  frequency  AQ. 
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V.     SINGLE  MODE  SOLUTIONS 

A.      SINGLE  MODE  MATRIX  SOLUTIONS 

In  the  case  of  a  spatially  complete  system  we  have  seen  that  the  SST  yielded 
frequency  independent  error  matrices  AK,  AM,  and  AC  that  could  be  used  to  correct  the 
stiflhess,  mass,  and  damping  matices  of  the  FE  model  in  such  a  manner  that  the  FRF  of  the 
corrected  FE  model  was  exactly  equal  to  that  of  the  experimental  system.  In  the  case  of  a 
spatially  incomplete  system  the  constituent  solutions  AK,  AM,  and  AC  given  in  chapters 
III  and  IV  were  in  general  frequency  dependent  solutions  which  serve  only  as  corrections 
to  the  reduced  FE  model.  Ideally  we  seek  frequency  independent  solutions  which  are 
corrections  to  the  full  FE  model.  For  now  we  shall  only  deal  with  the  simpler  problem  of 
trying  to  find  frequency  independent  solutions  which  serve  as  corrections  to  the  reduced 
FE  model. 

We  shall  employ  Equation  (2.30).  For  a  given  experimental  system  mode,  ©n, 
consider  a  frequency  bandwidth  [Qi,  Qu]  such  that  cOn€[Qi,  Qu].  Let  S={Qi  Q2,...,Om}  be 
a  frequency  sampling  of  the  bandwidth  [Qi,  Qu],  i.e.,  Qi<Q,<Qufor  each  QjeH.  For  each 
QieS  we  can  form  the  error  impedance  matrix  AZ(Oi).  We  apply  Equation  (2.30)  to  the 
partition  yeilding  the  following  system  of  m  equations  in  three  unknowns: 


AZC(Q,)" 

"/ 

-Qf/     j'Qj" 

~ak?-~ 

AZC(Q,) 

— 

/ 

-a]i   jnj 

aa/;- 

AC- 

AZc(Qj 

I 

-n2j  jnj 

(5.1) 


We  will  demonstrate  by  example  that  for  properly  chosen  bandwidths  [Qi,  Qu],  the 
solution  AK""t  AM,?" ,  AC""  of  Equation  (5.1)  approximately  corrects  the  Analytic  FRF 
over  the  frequency  bandwidth  [Qi,Qu],  i.e.,  \S.Ft(Q)  is  the  value  of  the  FRF  matrix  of  the 
experimental  system  at  a  frequency  Q,  //a(Q)the  value  of  the  FE  reduced  analytic  FRF 
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matrix  at  Q,   and  //"^(Q)  the  value  of  the  single  mode  corrected  reduced  Analytic  FRF 
matrix  at  Q  for  an  experimental  system  mode  (On  where  Qe  [Qi,Qu]  then 


H*(n)-HUn)^H*(n)-H°(n) 


(5.2) 


We  will  refer  to  the  solution  AK°n ,  AM®" ,  AC®"  as  a  single  mode  solution  at  Q)n. 

In  the  case  of  a  spatially  complete  system,  such  as  the  spatially  complete  beam  discussed  in 
chapter  III,  all  single  mode  solutions  are  found  to  be  identical  to  the  unique  frequency 
independent  solution  that  was  obtained  in  chapter  III,  hence  the  corrections  are  valid 
throughout  the  frequency  range  of  the  experimental  system  and  the  left  hand  side  of 
Equation  5.2  is  zero.  For  a  spatially  incomplete  system  this  is  not  the  case. 

In  support  of  Equation  (5.2)  we  chose  a  1  Hz  frequency  bandwidth  centered  on 
mode  1  (coi=25.21  Hz)  of  our  spatially  incomplete  beam  as  defined  in  Chapter  III.  For  our 
frequency  sampling,  E,  we  will  use  a  sampling  frequency  of  .5  Hz  which  results  in  3  points 
which  are  equally  spaced  over  the  interval,  E={24.7178  Hz,  25.2178  Hz,  25.7178  Hz}. 
Solving  Equation  (5.1)  we  get 


6K?~ 


-3.91      -13.60      21.70 
-13.60      170.5      -187.5 
-21.70    -187.5      194.6 


(5.3a) 


AAC'  = 


-0.0001      0.0006      -0.0005 

0.0006     -0.0033      0.0026 

-0.0005    -0.0026    -0.0016 


(5.3b) 


AC  = 


0  +  0.0270/  0-0.2291;  0  +  0.2277/ 
0-0.2291/  0  +  1.5078/  0-1.4627/ 
0  +  0.22777     0-1.4627;     0  +  1.3781/ 


(5.3c) 


Figure  5-1  shows  a  comparison  of  experimental,  uncorrected  and  mode  1 
corrected  FRFs  for  our  spatially  incomplete  beam  using  a  1  Hz  frequency  bandwidth 
centered  on  Mode  1  and  a  frequency  sampling  consisting  of  3  frequencies  equally  spaced 
over  the  bandwidth.  Figure  5-2  shows  the  results  of  including  mode  1  and  mode  2  in  the 
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frequency  bandwidth  [Q\,  Qu]-  As  Figure  5-2  shows  Equation  5.1  is  not  as  accurate  when 
multiple  frequencies  are  included  in  the  bandwidth.  Figure  5-3  is  a  comparison  plot  of  the 
experimental  and  analytic  FRF  versus  the  single  mode  matrix  solution  corrected  FRFs  over 
25,  10,  and  1  Hz  bandwidths  using  a  3  point  frequency  sampling.  Figure  5-3  shows  that 
the  precedure  is  sensitive  to  the  size  of  the  bandwidth  over  which  the  solution  is 
computed,  i.e.,  better  accuracy  is  achieved  with  smaller  bandwidths.  Figure  5-4  is  a 
comparison  plot  of  the  experimental  and  uncorrected  analytic  FRFs  versus  the  single  mode 
matrix  solution  corrected  FRFs  computed  over  a  1  Hz  bandwidth  using  frequency 
samplings  of  200,  50,  10,  and  3  points  equally  spaced  over  the  bandwidth.  Figure  5-4 
shows  that  the  procedure  is  fairly  insensitive  to  sampling  frequency. 

Figure  5-5  is  a  comparison  plot  of  the  experimental  and  analytic  FRFs  versus  the 
single  mode  solution  corrected  FRF  for  a  1  Hz  bandwidth  using  a  3  point  frequency 
sampling  in  the  case  of  a  spatially  complete  beam.  As  Figure  5-5  shows  the  procedure  is 
exact  for  spatially  complete  systems. 


46 


1-2 

»— I 
60 
O 

H 

1-4 


x 


-8 


0 


if 

If 

>..!» i 

1 

s ;..... 

I 

T  ' 

_ 

mental  FRF 
1  Corrected  FRF 

IP-  - 

bxpcn 
Mode 

.  " If  i  .  .  .  A 

'                I                ' 

i 

A!1     AS ;  'N\  ■ 

.'....;! .... 

<^r, 

V 

i 

■  - 

1 

1 

1 

i 

200 


400 


600 


800 


1000 


1200 


10 


I  - 

60 

o 

I 


-10 


1 

1       ■■ ■■  ■ 

i                 i  

i 

i 

>  ■  —                                    ^>- 

|/V" : 

V 

1 

i                  i                  i 

• 

10 


15 


20 


25  30 

Omega  (Hz) 


35 


40 


45 


50 


Figure  5-  1    Experimental  and  uncorrected  Analytic  FRFs  vs  single  mode  solution  at 
mode  1  corrected  FRF  using  a  3  point  frequency  sampling  of  a  1  Hz  bandwidth. 
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Figure  5-  2   Experimental  and  uncorrected  analytic  FRFs  vs  single  mode  corrected  FRF 
using  a  15  point  frequency  sampling  of  a  bandwidth  that  includes  modes  1  and  2. 
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Figure  5-  3    Spatially  incomplete  experimental  FRF  vs  single  mode  matrix  solutions  at 
mode  1  using  3  point  frequency  samplings  of  25,  10,  and  1  Hz  bandwidths. 
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Figure  5-  4    Spatially  incomplete  experimental  FRF  vs  single  mode  matrix  solutions  at 
mode  1  using  3,  10,  50,  and  200  point  frequency  samplings  of  a  1  Hz  bandwidth. 
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Figure  5-  5    Spatially  complete  experimental  FRF  vs  single  mode  matrix  solutions  at 
mode  1  using  a  3  point  frequency  samplings  of  a  1  Hz  bandwidth. 
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B.     SINGLE  MODE  INTEGRAL  SOLUTIONS 

In  what  follows  we  wish  to  employ  the  integral  formulas  of  reference  (5)  to  obtain 
a  results  similar  to  that  of  Equation  (5.2).  To  accomplish  this  we  shall  need  to  recast  the 
impedance  equations  of  Chapter  II  in  terms  of  velocity.  We  start  with  the  equation  of 
motion  of  a  general  2nd  order  linear  system 

Kx  +  Cx  +  Mx  =  f  (5.4a) 

where 

x  =  Xejnt  (5.4b) 

f=FejQl  (5.4c) 

If  we  differentiate  Equation  (5.4b)  we  get  that 

x  =  jOXe**  =  jQx  (5.5) 

hence 

x  =  -£r  (5-6) 

jQ, 

differentiating  Equation  (5.4b)  twice  we  get  that 

x  =  jdyQXe^)  =  jQx  (5.7) 

substituting  Equations  (5.6  )and  (5.7)  into  Equation  (5.4a)  we  get 

K—  +  Cx  +  MjQx  =  f  (5.8) 

jQ 

We  can  rewrite  equation  (5.8)  as 

(  K  \ 

—  +  C  +  Mjnx  =  f  (5.9) 

\jQ  J 

Following  reference  (5)  we  can  take  the  Laplace  transform  of  equation  (5.9)  and 
thus  write  the  impedance  of  the  linear  system  in  terms  of  the  complex  Laplace  parameter  s 
as 
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Z(s)  =  Ms  +  C  +  —  (5-10) 

s 

In  reference  (5)  the  impedance  matrix  of  a  general  system  is  represented  as  an 
infinite  Laurent  series  about  the  orgin  in  powers  of  the  complex  Laplace  parameter  s  as 

Z(s)=  IV  =  Limn_SAnsn+-+Als  +  A,  +  A_±+-+A_n\)  (5.11) 

n=-oo  S  S 

After  defining  a  truncated  Laurent  series 

Z(s)  =  Als  +  A0+A  l-  =  Ms  +  C  +  —  (5.12) 

s  s 

which  approximates  the  impedance  matrix,  an  error  function  E(s),  and  a  cost  J  are  defined 
as 

E(s)  =  Z{s)-Z(s)  (5.13) 

and 

J=iJE(s)\\E\W(s)-*\  (5-14) 

where  W(s)  is  a  complex  valued  weighting  function,  the  subscript  E  denotes  the  euclidian 
norm  of  a  matrix  and  the  integration  is  performed  over  a  prescribed  path  P  in  the  complex 
plane.  By  setting  the  partial  derivatives  of  J  with  respect  to  M ,  K ,  and  C  to  zero  the 
authors  of  reference  (5)  obtained  expressions  for  matrices  M ,  K ,  andC  which 
approximate  the  stiffness,  mass  and  damping  matrices  about  a  mode  of  the  system.  The 
expression  for     M  ,K ,  and  C  are  as  follows: 


M=      l 


ab-c' 


>  1 

b  JQZj(iQ)\W(iQ)\dQ  -  c  J— Z7(iQ)|ir(iQ)|dQ 

(5-15a) 
=  -  JZa(/Q)|JF(/ft)|dQ  (5-15b) 


C=- 
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r-    l 


ab-c' 


n2  o> 

c  JQZj(iQ)\W(iQ)\dQ  -  a  J— Z7(iQ)|JF(/Q)|dn 

(5- 15c) 


where 


"2 

a=  JQ2\W(iQ)^Q  (5-15d) 


«2       1 

/3=  J— |Fr(/Q)|dQ  (5-15e) 

n2 
c=  J|JT(/Q)|dn  (5-15f) 

Z(/Q)  =  Z^/Q)  +  /Z7(/Q)  (5-15g) 

We  will  apply  Equation  (5.15)  to  AZ  as  defined  by  the  SST  and  obtain 
matrices  AK ,  AM ,  and  AC  which  will  serve  as  correction  matrices  of  the  FRF  of  the 
system  about  a  given  mode  of  the  system. 

As  an  example,  let  us  use  Equation  (5.15)  to  compute  the  (1,1)  element  of  the 
correction  matrix  AC  at  mode  1  (coi=25.21  Hz)  for  our  spatially  incomplete  system.  We 
take  the  weighting  function  to  be  W(iQ)=l/(iO)  ,  the  path  P  is  taken  to  be  along  the 
imaginary  axis  from  155.30  rads/sec  to  161.58  rads/sec. 

To  compute  this  line  integral  we  shall  use  the  trapezoid  rule  with  a  three  point 
frequency  sampling 

S  =  {\5530rads  I  sec,158.44rafr  /  sec,161.58rafc  /  sec}  (5.16) 

of  the  straight  line  path  P  (sampling  frequency  of  n  rads/sec).  We  first  compute  the  value 
of  the  constant  a  of  Equation  (5.15d).  Using  vector  notation  we  have  that 

Q(S)  =  [\553Qrads  I  sec,158.44ra<fr  /  sec,161.58/rafr  /  sec]  (5.17a) 

and 

JV{E)  =  [0  -  0.0064/ sec/  rads,0  -  0.0063y  sec/  rads,0  -  0.0062/ sec/  rads]      (5. 1 7b) 
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hence  the  integrand  is  the  vector 

Q2(E)-|^(E)|  =  [155.30,158.44,161.58]  (5.18) 

Using  the  trapezoid  rule  to  compute  the  integral,  we  have  that 

a=  7!n2(E)|(r(S)|^  =  ^(^  + 158.44 +  i^il  =995.S6"*y  ,  (5.19) 

155.30  L  Z  Z 

b=. 00000158  sec/rads2  and  c=0.0397  are  computed  in  a  similar  manner. 

To  actually  compute  the  matrix  AC  ,we  use  the  SST  to  compute  the  matrices 
AZ(l 55.30  rads/sec),  AZ(1 55.44  rads/sec),  and  AZ(161.58  rads/sec)  for  the  points  of  the 
frequency  sampling  of  the  straight  line  path  P.  Using  Equation  (2.26)  the  matrix 
AZ(155.30  rads/sec)  is  seen  to  be 


AZ(l55.44^%c)  = 


"(0.64£  -  9)  +  0.02 j    (-0.04£  -  8)  -  0.04/     (0.03£  -  7)  +  0.01; 
(-0AE  -  8)  -  0.04./     (-0.4£  -  9)  -  0.04/      (0.25£  -  7)  -  0.1 0/ 
(0.3£  -  8)  +  0.0 1/       (0.25i:  -  7)  -  0.1 0/     (-0.24E  -  7)  +  0.15;_ 

(5.20) 

The  integral  in  Equation  (5.15b)  is  computed  on  an  element  by  element  basis.  For 

the  (1,1)  element  of  Equation  (5.15b),  if  we  collect  the  (1,1)  elements  of  the  matrices 

AZ(1 55.30  rads/sec),    AZ(1 55.44   rads/sec),    and   AZ(161.58   rads/sec)    and   use   the 

weighting  vector  given  in  Equation  (5.17b)  the  integrand  of  Equation  (5.15b)  is  the  vector 

AZiU)(S)|FT(S)|  =  [0.0064    0.040    -0.848][0.64    0.63    0.62](£  -  9)  lb%n  se(/rads 

(5.21a) 
AZiU)(S)|^(E)|  =  [0.0041     0.2885    0.0036](£-  ^)lhf/in^/rads  (5.21b) 

Computing  the  integral  we  have 

(5.22a) 
AC(1,1)  =  0.0232£  -  6lbf/iri  (5.22b) 

Performing  the  above  procedures  for  the  remaining  elements  of  Equation  (5.15b)  we  get 
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/in 


AC  = 


'  0.0232     -0.0443     0.0300 " 

-0.0443      0.1278      -0.0914 

0.0300     -0.0914      0.0645 


M;% 


(5.23a) 


In  a  similar  manner  we  have 


AM  = 


-0.0001  -0.0001  0.0002 ' 
-0.0001  0.0015  -0.0020 
0.0002     -0.0020     0.0027 


Ibm 


AK  = 


-6.0476      4.5347        3.6757 
4.5347      51.2196     -71.7061 
3.6757     -71.7061     85.4417 


Ibfj 


in 


(5.23b) 


(5.23c) 


Figure  5-6  is  a  comparison  plot  of  experimental  and  uncorrected  analytic  FRFs 
versus  single  mode  integral  solutions  at  mode  1  using  a  frequency  sampling  of  3  points 
equally  spaced  over  bandwidths  of  25,  10,  and  1  Hz.  As  with  the  matrix  solutions,  the 
integral  solutions  are  sensitive  to  the  length  of  the  bandwidth  over  which  the  solutions  are 
computed  with  smaller  bandwidths  yielding  more  accurate  solutions.  Figure  5-7  is  a 
comparison  plot  of  experimental  and  uncorrected  analytic  FRFs  versus  single  mode 
integral  solutions  at  mode  1  over  a  1  Hz  bandwidth  using  frequency  samplings  of  3,10,50, 
and  200  points  from  the  bandwidth.  As  with  the  matrix  solutions,  the  integral  solutions  are 
also  insensitive  to  the  sampling  frequency. 

Figures  5-8,  5-9,  and  5-10  are  comparison  plots  of  matrix  and  integral  solutions  at 
mode  1  computed  over  25,  10,  and  1  Hz  bandwidths  respectively,  using  three  point 
frequency  samplings.  These  three  plots  show  that  for  these  frequency  sampling,  the  matrix 
solutions  are  more  accurate  than  the  integral  solutions.  We  must  note  that  we  have  used 
the  trapezoid  rule  for  computation  of  the  integrals.  It  is  expected  that  better  accuracy 
would  be  achieved  from  the  integral  solutions  if  a  higher  order  integration  method  such  as 
Simpson's  rule  were  used. 
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Figure  5-  6    Spatially  incomplete  experimental  FRF  vs  single  mode  integral  solutions  at 
mode  1  using  3  point  frequency  samplings  of  25,  10,  and  1  Hz  bandwidths. 
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Figure  5-7    Spatially  incomplete  experimental  FRF  vs   single  mode  integral  solutions  at 
mode  1  using  3,  10,  50  ,  and  200  point  frequency  samplings  of  a  1  Hz  bandwidth. 
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Figure  5-  8   Experimental  FRF  vs  single  mode  integral  and  matrix  solutions  at  mode  1 
using  a  25  Hz  bandwidth  with  a  3  point  frequency  sampling. 
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Figure  5-  9   Experimental  FRF  vs  single  mode  integral  and  matrix  solutions  at  mode  1 
using  a  10  Hz  bandwidth  with  a  3  point  frequency  sampling. 
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Figure  5-  10   Experimental  FRF  vs  single  mode  integral  and  matrix  solutions  at  mode  1 
using  a  1  Hz  bandwidth  with  a  3  point  frequency  sampling. 
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VI.     MULTIPLE  MODE  SOLUTIONS 

A.      MULTIPLE  MODE  MATRIX  SOLUTIONS 

We  wish  to  extend  the  results  of  Chapter  V  to  multiple  modes  of  the  system  under 
consideration.  To  accomplish  this  we  shall  simply  extend  Equation  (5.1)  over  multiple 
frequency  samplings.  For  a  set  of  system  modes  {©i,  ci>2,...,c0n},  let  H(cOi)={Qii, 
Qi2,...,Qik}  be  a  frequency  sampling  of  a  bandwidth  [Qii,QjU]  about  co,  for  i=l,2,...,n.  We 
shall  apply  Equation  (2.30)  to  the  frequency  sampling 


B-UW-.) 


(6.1) 


which  is  the  set  theoretic  union  of  the  sets  H(coj)  for  i=l,2,...,n.  If,  for  simplicy,  we  restrict 
ourselves  to  equally  sized  frequency  samplings  we  obtain  a  set  on  nk  equation  in  three 
unknowns 

7    -Q2UI    jQnf 


AZc(Qn) 

azc(q,) 

AZC(Q„,) 

One  can  solve  Equation  (6.2)  for 


/    -Q2nkI    jQ.nkI 


am: 


AC; 


(6.2) 


am; 

AC? 


as  we  have  done  in  Chapter  V  using  the 


familiar  MATLAB  puedoinverse  function  but  better  results  are  obtained  if  we  weight  the 
equations  as  discussed  in  [Ref.  6].  A  good  weighting  is  to  assign  the  weight  co*  to  those 
equations  associated  with  the  frequency  sampling  H(o)i).  Alternately  we  could  assign  the 
weight  l/coi  to  those  equations  associated  with  the  frequency  sampling  E(coj).  The  ©i 
weighting  results  in  the  solution  being  more  accurate  at  the  lower  modes  while  the  1/cd; 
weighting  gives  better  accuracy  at  the  higher  modes. 

We  note  that  the  solution  to  Equation  (6.2)  approximately  corrects  the  analytic 
FRF  in  the  sense  of  Equation  (5.1)  when  the  matrices  AZ(Dtj)  as  given  by  Equation  (2.26), 
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are  computed  over  a  'good'  error  set,  Cen-  In  general,  due  to  the  smearing  that  occurs 
during  reduction  of  the  full  analytic  system  to  the  reduced  analytic  system,  the  set  Cen  is 
not  the  same  set  as  the  intersection  of  the  spatially  complete  error  set  and  the  A  set.  It  is  in 
general  a  larger  set  than  this  intersection  and  the  size  is  a  function  of  the  type  of  error,  i.e., 
mass,  damping,  or  stiffness  error  and  the  locations  of  the  errors.  Figure  6-1  is  a 
comparison  plot  of  experimental  and  uncorrected  analytic  FRFs  versus  the  multiple  mode 
matrix  solution  corrected  FRF  with  the  solution  having  been  computed  over  modes  1  and 
2  using  a  1  Hz  bandwidth  about  each  mode  with  a  3  point  frequency  sampling  over  the 
bandwidth  at  each  of  the  modes.  Figure  6-2  is  a  comparison  plot  using  a  matrix  solution 
computed  over  mode  1  through  4  again  using  a  1  Hz  bandwidth  about  each  mode  with  a  3 
point  frequency  sampling  over  the  bandwidth  at  each  of  the  modes. 
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Figure  6- 1  Experimental  FRF  vs  multiple  mode  matrix  solutions  at  modes  1  and  2  using 
a  1  Hz  bandwidth  with  a  3  point  frequency  samplings  at  each  mode. 
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Figure  6-  2   Experimental  FRF  vs  multiple  mode  matrix  solutions  at  modes  1  through  4 
using  a  1  Hz  bandwidth  with  a  3  point  frequency  samplings  at  each  mode. 
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B.     MULTIPLE  MODE  INTEGRAL  SOLUTIONS 

We  wish  to  extend  single  mode  integral  solutions  to  multiple  modes  of  the  system. 
As  with  the  multiple  mode  matrix  solutions  we  simply  take  as  the  path  of  integration  a 
path  which  is  the  set  theoretic  union  of  suitable  paths  at  each  of  the  modes  under 
consideration.  As  with  the  matrix  solutions  one  can  chose  the  weighting  function  W(s)  to 
be  1/(Q)  or  Q  to  achieve  lower  or  higher  mode  accuracy.  As  stated  in  the  previous  section 
it  is  required  that  a  'good'  error  set  is  known.  Figure  6-3  is  a  comparison  of  experimental 
and  uncorrected  FRF  versus  multiple  mode  integral  solution  corrected  FRF  with  the 
solution  having  been  computed  over  modes  1  and  2  using  a  1  Hz  bandwidth  at  each  mode 
with  a  3  point  frequency  sampling  over  the  bandwidth  at  each  of  the  modes.  Figure  6-4  is 
a  comparison  using  a  solution  computed  over  mode  1  through  4  again  using  a  1  Hz 
bandwidth  at  each  mode  with  a  3  point  frequency  sampling  over  the  bandwidth  at  each  of 
the  modes.  Figure  6-5  is  a  comparison  plot  of  the  multiple  mode  matrix  and  integral 
solutions  over  modes  1  through  mode  4  using  a  1  Hz  bandwidth  at  each  mode  with  a  3 
point  frequency  sampling  over  the  bandwidth  at  each  of  the  modes. 
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Figure  6-  3  Experimental  FRF  vs  multiple  mode  integral  solutions  at  modes  1  and  2  using 
a  1  Hz  bandwidth  with  a  3  point  frequency  samplings  at  each  mode. 
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Figure  6-  4  Experimental  FRF  vs  multiple  mode  integral  solution  at  modes  1  through  4 
using  a  1  Hz  bandwidth  with  a  3  point  frequency  samplings  at  each  mode. 
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Figure  6-  5  Experimental  FRF  vs  multiple  mode  matrix  and  integral  solutions  at  modes  1 
and  2  using  a  1  Hz  bandwidth  with  a  3  point  frequency  samplings  at  each  mode. 
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C.     SINGLE  POINT  MULTIPLE  MODE  SOLUTIONS 

The  results  of  sections  A  and  B  can  be  performed  using  partitions  consisting  of  a 
single  point.  Figure  6-6  shows  the  results  of  performing  a  multiple  mode  matrix  solution 
over  the  first  three  modes  of  our  spatially  incomplete  beam  using  a  single  point  partition  at 
each  of  the  modes  under  consideration.  Figure  6-7  shows  the  results  of  performing  a 
multiple  point  integral  solution  over  the  first  three  modes  of  our  spatially  incomplete  beam 
using  a  single  point  partition  at  each  of  the  modes  under  consideration.  Figure  6-8 
compares  the  multiple  mode  matrix  and  integral  solution  computed  over  the  first  three 
modes  of  our  spatially  incomplete  beam  with  the  analytic  FRF.  In  figure  6-9  we  compare  a 
multiple  mode  matrix  solution  computed  over  the  first  three  modes  of  a  spatially  complete 
beam  using  a  single  point  partition  at  each  of  the  modes  with  the  beam's  spatially  complete 
analytic  FRF.  In  figure  6-10  we  compare  a  multiple  mode  integral  solution  computed  over 
the  first  three  modes  of  a  spatially  complete  beam  using  a  single  point  partition  at  each  of 
the  modes  with  the  beam's  spatially  complete  analytic  FRF.  Figure  6-11  compares  multiple 
mode  matrix  and  intgral  solution  computed  over  the  first  three  modes  of  a  spatially 
complete  beam  with  the  beam's  spatially  complete  analytic  FRF.  Figures  6-9  and  6-10 
show  that  .the  single  point  multiple  mode  solutions  are  exact. 
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Figure  6-  6  Experimental  FRF  vs  multiple  mode  matrix  solution  computed  at  modes  1 
through  3  using  a  1  Hz  bandwidth  with  a  single  point  frequency  samplings  at  each  mode 
for  a  spatially  incomplete  beam. 
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Figure  6-  7  Experimental  FRF  vs  multiple  mode  integral  solution  computed  at  modes  1 
through  3  using  a  1  Hz  bandwidth  with  a  single  point  frequency  samplings  at  each  mode 
for  a  spatially  incomplete  beam. 
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Figure  6-  8  Experimental  FRF  vs  multiple  mode  matrix  and  integral  solutions  computed  at 
modes  1  through  3  using  a  1  Hz  bandwidth  with  single  point  frequency  samplings  at  each 
mode  for  a  spatially  incomplete  beam. 
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Figure  6-  9  Experimental  FRF  vs  multiple  mode  matrix  solution  computed  at  modes  1 
through  3  using  a  1  Hz  bandwidth  with  a  single  point  frequency  samplings  at  each  mode 
for  a  spatially  complete  beam. 
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Figure  6-  10  Experimental  FRF  vs  multiple  mode  integral  solution  computed  at  modes  1 
through  3  using  a  1  Hz  bandwidth  with  a  single  point  frequency  samplings  at  each  mode 
for  a  spatially  complete  beam. 
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Figure  6-11  Experimental  FRF  vs  multiple  mode  matrix  and  integral  solution  computed 
at  modes  1  through  3  using  a  1  Hz  bandwidth  with  a  single  point  frequency  samplings  at 
each  mode  for  a  spatially  complete  beam. 
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300 


Figure  6- 12  Experimental  FRF  vs  multiple  mode  matrix  solutions  computed  over  1,  2,  3, 
and  4  modes  using  a  1  Hz  bandwidth  with  single  point  frequency  samplings  at  each  mode 
for  a  spatially  complete  beam. 
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Vn.     CONCLUSIONS  /  RECOMMENDATIONS 

A.  SUMMARY 

Frequency  Domain  Structural  Identification  using  the  structural  synthesis 
transformation  was  performed  on  a  simulated  experimental  free-free  beam  model.  The 
identification  was  performed  for  both  the  spatially  complete  and  the  spatially  incomplete 
case. 

•  SST  based  structural  identification  provides  an  exact  solution  for  the 
identification  of  FE  modeling  errors  given  spatially  complete  data. 

•  For  spatially  incomplete  systems  SST  based  structural  identification  provides 
a  frequency  dependent  solution  which  is  not  suitable  for  finite  element  modeling  error 
correcting. 

•  Using  both  matrix  and  integral  based  techniques  the  SST  can  provide  single 
mode  solutions  which  are  frequency  independent  correction  matrices  AK,  AM,  AC 
which  approximately  corrects  the  reduced  FE  model  within  a  frequency  bandwidth  of 
the  experimental  system  mode  under  consideration. 

B.  CONCLUSIONS 

This  thesis  has  clearly  demonstrated  that  for  spatially  incomplete  systems,  single 
mode  frequency  independent  solutions  can  be  found  which  correct  the  reduced  finite 
element  model  in  a  neighborhood  of  a  given  mode  of  the  experimental  system.  It  has  also 
been  shown  that  the  concept  of  a  single  mode  solution  can  be  expanded  to  that  of  multiple 
mode  solutions  which  are  frequency  independent  correction  matrices  which  approximately 
corrects  the  reduced  FE  model  throughout  a  frequency  bandwidth  which  includes  more 
than  1  mode  of  the  experimental  system. 

C.  RECOMMENDATION 

The  purpose  of  this  thesis  was  to  find  frequency  independent  multiple  mode 
solutions  which  could  be  used  to  approximately  correct  reduced  finite  element  models. 
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Although  satisfactory  results  were  obtained,  investigation  is  still  required  in  the  following 
areas: 

•Determine  the  relationship  between  errors  in  a  full  FE  model  and  those  of  an 
associated  reduced  FE  model. 

•Determine  a  method  to  'pullback'  reduced  FE  model  correction  to  full  order 
FE  model  corrections. 

•Investigate  use  of  the  integral  formulation  of  reference  (5)  in  analysis  of  the  O- 
set  system. 
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APPENDIX 

The  following  is  a  brief  description  of  MATLAB  routines  employed  in  this  thesis: 

•SST.M   -   Generates   MAT   files   containing   spatially   complete   or   spatially   incomplete 
experimental  FRF,  analytic  FRF,  localization  matrix,  and  SST  solutions. 

•PLOTSST.M  -  Uses  files  produced  by  SST.M  to  generate  plots  used  in  chapters  3  and  4. 

•CHAP3.M  -  Sets  parameters  and  calls  SST.M  and  PLOTSST.M  to  generate  figures  for 
chapter  3 . 

•CHAP4.M  -  Sets  parameters  and  calls  SST.M  and  PLOTSST.M  to  generate  figures  for 
chapter  4. 

•CHAP5.M  -  Generate  figures  for  chapter  5. 

•CHAP6.M  -  Generate  figures  for  chapter  6. 

•SETUP.M  -  Generates  MAT  files  for  FE  models. 

•BEAMMDL.M  -  Setup  FE  models. 

•FSTATIC.M  -  Perform  static  reduction  of  mass  and  stiffness  matrices. 

•FIRS_TAM.M  -  Perform  IRS  reduction  of  mass  and  stiffness  matrices. 

•FREQMODE.M  -  Returns  FE  model  frequencies. 

•NDX3D.M  -  Indexes  a  series  of  2-D  matrices  into  a  single  3-D  matrix. 

•INTSUB.M  -  Decomposes  impedance  matrix  into  mass,  stiffness,  and  damping  matrices  using 
integral  techniques. 

•MYTRAP.M  -  Computes  integrals  using  trapezoidal  method 
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%%%%%%%%%%%%0/o0/o%0/o0/o%0/o%%%%%%0/o%%%%%%%%0/o%%%0/o 

%%%%%Generalsetup%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%0/o%0/o%%0/o%0/o%%%0/o0/o%%%%%%0/o%0/o0/o%0/o% 

clear,clc         ;%%clear  workspace. 
5        closeall  ;%close  any  open  figure  windows 

j=sqrt(-l)        ; 

load  sstconf     ;%load  A-set  and  O-set 

aset(oset)=[]; 
1 0        aset_size=length(aset); 

oset_size=length(oset); 

load  beamdata 

if  length(oset)==  0 
complete=l; 
1 5        else 

complete=0; 

end 

%%semppl0tlabels%%%%%%%%%%%%%%%%%%%%%%%%%%%0/0%%%%%%%%%%% 
20        if  (oset==[])&(posofMasserr(  1  )~^0)&(posofStifferT(  1  >^0) 
casenaine=('Spatially  Complete  Mass  &  Stiffness  error'); 
elseif  (oset=[])&(posofMasserr(  1  )~=0) 

casename=('Spatially  Complete  Mass  error1); 
elseif  (oset==[])&(posofStifferr(l  y=0) 
25  casename=('Spatially  Complete  Stiffness  error'); 

elseif  (posofMasserr(  1  )~^0)&(posofStifferr(  1  )~=0) 

casename=('Spatially  Incomplete  Mass  &  Stiffness  error'); 
elseif  (posofMasserr(  1  )~=0) 
casename=('Spatially  Incomplete  Mass  error1); 
30        elseif  (posofStifrerr(l)-^O) 

casename=('Spatially  Incomplete  Stiffness  error"); 
end 


35 


%%get  reduced  A  set  &  O  set  frequencies%%%%%%%%%%%%% 


kexto=k_anal(oset,oset);  %%  stiffness 
mexto=m_anal(oset,oset);  %%  and  mass  matrices 
if  complete 
kstat=k_anal; 
40  mstat=m_anal; 

else 
if  static  ==  0 

[kstat,mstat]=fstatic(k_anal,m_anal,oset,aset);%%  get  reduced  K  &  M 
elseif  static  ==  1 
45  [kstat,mstat]=firs_tam(k_anal,m_anal,oseLaset);%%  get  reduced  K  &  M 

else 

kstat=k_anal(aset,aset); 
mstat=m_anal(aset,aset); 
end 
50        end 

cstat=sqrt(- 1  )*struc_damping.  *kstat; 
[u,lambdaa,c]=freqmode(k_anal,m_anal)     ;%get  freqs 
[u,lambdared,c]=fxeqmode(kstat^nstat) ;% 
[u,lambdax,c]=freqmode(k_exp,m_exp)       ;% 
55         [u,lambdaexto,c]=freqmode(kexto,mexto);% 
omegaa=sqrt(  lambdaa); 
omegax=sqrt(lambdax); 
omegared=sqrt(lambdared); 
omegaexto=sqrt(lambdaextoX 
60        for  i=  1 .4     %  tins  to  delete  lixed  bod}  modes 
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if  real(omegax(l))  <  10A(-3) 

omegax=oinegax(2:length(omegax)); 
end 

if  reaI(omegaa(l))  <  10A(-3) 
65  omegaa=omegaa(2:length(omegaa)); 

end 
ifcomplete==0 

ifreal(omegared(l))<  10A(-3) 
omegared=omegared(2:length(omegared)); 
70  end 

ifreal(omegaexto(l))<  10A(-3) 

omegaexlo=omegaexlo(2:length(omegaexlo)); 
end 
end 
75        end 

%%%%%fmd  true  mass  and  stiffiiess  errors%%%%%%%%%%%%%%%%%%% 
lrue_stiffriess=k_exp-k_anal; 
stiffiiess_cset=find(diag(true_stiffhess)); 
80        true_iuass=m_exp-m_anal; 

mass_cset=find(diag(true_mass)); 
true_damping=c_exp-c_anal; 
dainping_cset=find(diag(true_damping)); 
if  length(mass_cset)  >  0 
85  for  i=  1  :length(mass_cset) 

x=find(stiffhess_cset==mass_cset(i)); 
if~<length(x)>0) 

stiffness_cset=[stiffness_cset;  mass_cset(i)]; 
end 
90  end 

end 

if  length(damping_cset)  >  0 
for  i=l:length(damping_cset) 

x=fmd(sUffness_cset=~damping_cset(i)); 
95  if~<length(x)>0) 

stiffiiess_cset=[sliffness_cset;damping_csel(i)]; 
end 
end 
end 
100        true_cset=sort(stiffness_cset) 

save  true_errs  true_cset  true_mass  true_sliffness  true_damping 
clear  truemass  true_stiffhess  truedamping 

clear  area  eeii  force  g  gc  gcx  gk  gkx  gin  ginx  goblc  goblcx  goblk  goblkx 
clear  goblm  goblmx  ke  lainbdaa  lambdaexto  lainbdared  lambdax  luinpdainp 
1 05        clear  lumpmass  lunipspring  rae  pho  posofDaniperr  posofMasserr  posotEtifferr 
clear  stifihess_cset 

%%%setup  frequency  raige  for  midysis%%%%%%%%%%%%%%%%%%%%%%0/o0/o0/o0/o%0/o0/o%0/o0/o 
freqtop=omegax(luglunode>i-.  1  *omegax(luglunode); 
1 1 0        freqbottom=omegax(lowmode)-.  1  *omegax(lowmode); 
freqbottom=10*2*pi;    %start  at  1 0  hz 
freqtop=  1 200*2  *pi;     %end  at  1 200  hz 
%%number  of  points  to  plot 

1 1 5        nuinpomts=400;%%%%%SET  NUMBER  OF  POINTS  TO  USE  FOR  THE  SIMULATION 
fineness=numpoints; 

vv=freqboltom:(freqtop-freqbottom)/(nunipoints-l):freqtop; 

I20        %%LorateCset%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%0/o%0/o%0/o0/o%0/o%0/o 
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keepgoing-n' 

wfreq=omegax(lowmode) ; 
nextinode=lowmode; 
tol=l  ; 

125        defaultstepsize=.l; 

while  (keepgoing='n') 
z_anal_red=kstat+j  *wfreq*cslat-wfreqA2  *instat ; 
h_anal_red=inv(z_anal_red)  ; 

z_exp=k_exp+j  *wfreq*c_exp-wfreqA2  *m_exp    ; 
1 30  h_exp=inv(z_exp)  ; 

h_exp=h_exp(aset,aset) 

L=z_anal_red*(h_anal_red-h_exp)*z_anal_red  ; 

ci=fmd(abs(diag(L))>tol)  ; 

lemp_cset_size=lengtli(ci); 
135  iflength(ci)>=l 

cset_size=length(ci)  ; 

figure(l) 

plot(aset,  abs(diag(L)))  ; 

hold  on 
1 40  plot(aset(ci),zeros(  1  ,length(ci)),'x') 

tallest=max(abs(diag(L)))       ; 

title([TOL  =  *,nuiii2str(tol),'  Cset  size  =  ',int2slr(cset_SLze),'  Omega  =  ',num2str(wfreq/(2*pi)),'  Hz']) 
ylabel('lbf/in') 
if  complete 
145  xlabel(T)OF) 

else 

xlabel('ASET  DOF) 
end 

hold  off 
150  k=menu(  '  Choose  an  action  ',.•• 

'hicrease  Tolerance  ',.•• 

'Decrease  Tolerance ',... 
Increase  frequency ',... 
Decrease  frequency ',... 
155  '    Previous  mode    ',... 

'    Next  Mode 
'Set  Print  switch  on',... 
'  hicrease  stepsize ',... 
'  Decrease  stepsize ',... 
160  'Go  '); 

if  k==l 
stepsize=defaultstepsize; 
ctr=0; 

while  length(ci)  >=  temp_cset_size 
165  ctr=ctr+l; 

ifctr==10 
stepsize=stepsize*  1 0; 
ctr=0; 
end 
170  tol=min(abs(tol+stepsize*(tol)),tallest) ; 

ci=fmd(abs(diag(L))>lol) 
iflengtli(ci)>=  1 
cset_size=length(ci)  ; 

figure(l) 
175  plot(aset,  abs(diag(L)))  ; 

hold  on 

plot(aset(ci),zeros(  1  ,length(ci)),'x') 
tallest=max(abs(diag(L)))       ; 

title([TOL  =  ',num2slr(tol)/  Csel  size  =  '.int2str(cse(_size),'  Omega  =  ',num2str(ufreq/(2*pi)).'  Hz']) 
180  ylabek'lbf/m') 
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if  complete 

xlabel('DOF) 
else 

xlabel('ASET  DOF) 
185  end 

hold  off 
end 
end 
elseif  k==2 
1 90  stepsize=defaullstepsize; 

ctr=0; 

while  length(ci)  <=  temp_cset_size 
ctr=clrH; 
ifctr==10; 
195  stepsize=stepsize  *2 ; 

cti=0; 
end 
tol=max(abs(tol-stepsize*(tol)),.00001); 
ci=fmd(abs(diag(L))>tol) 
200  iflength(ci)>=  1 

cset_size=lengdi(ci)  ; 

figure(l) 

plot(aset,  abs(diag(L)))         ; 
hold  on 
205  plot(aset(ci),zeros(  1  ,length(ci)),'x') 

lallest=max(abs(diag(L)))      ; 

tille([TOL  =  ',num2str(tol),'  Cset  size  =  ',int2str(cset_size),'  Omega  =  ',num2str(wtreq/(2*pi)),' 
Hz']) 

ylabel('lbf/in') 
210  if  complete 

xlabel(DOF) 
else 

xlabel('ASET  DOF') 
end 
215  hold  off 

end 
end 
elseif  k==3 
while  length(ci)  ==  temp_cset_size 
220  wfreq=mm(wfreq+.00001  *(freqtop-wfreq),freqtop); 

z_anal_red=kstat+j*w(count)*cstat-wlreqA2*mstat ; 
h_anal_red=inv(z_anal_red)  ; 

z_exp=k_exp+j  * wfreq*c_exp-\vfreqA2  *m_exp    ; 
h_exp=inv(z_exp) 
225  h_exp=h_exp(aset,aset)  ; 

L=z_anal_red*(h_anal_red-h_exp)*z_anal_red  ; 

tallest=max(abs(diag(L)))       ; 
tol=min(tol+.001  *tallest,lallest) ; 
ci=fmd(abs(diag(L))>tol)  ; 

230  cset_size=lengdi(ci)  ; 

figure(l) 

plot(aset,  abs(diag(L)))  ; 

hold  on 

plot(aset(ci),zeros(  1  ,lengtli(ci)),'x') 
235  title([TOL  =  ',nmn2str(tol),'  Cset  size  =  ',mt2str(cset_size)/  Omega  =  ',num2str(wfreq/(2*pi)),'  Hz']) 

ylabel('lbf/in') 
if  complete 

xlabel('DOF') 
else 
240  xlabel('ASET  DOF) 
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end 

hold  off 
end 

elseif  k==4 
245  wfreq=max(wfreq-.05*(wlreq-freqbottom),freqboUom); 

elseif  k==5 

nextniode=max(lowruode,nextmode-l ); 
wfreq=omegax(nextinode); 
elseif  k==6 
250  nexunode=min(lughmode,nextmode+ 1 ); 

wfreq=omegax(nextmode); 
elseif  k==7 

pswitch=y  ; 

elseif  k==8 
255  defaultstepsize=defaultstepsize/10; 

elseif  k==9 

defaultstepsize=defaullstepsize*  1 0; 
else 

keepgoing=y  ; 

260  cset=aset(ci); 

cset_rel=ci  ; 

end 
else 

tol=tol*.9; 
265  end 

end 

close(l) 

save  extrtset  cset  cset_rel 
270        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%0/o%0/o 
%%% 

%%%%%%%%%%%%%%%%%%FORCETHECSET%%%%%%%%%%%%%%%%%%%%%%%%%0/o%0/o0/o0/o% 
%cset=[7  9  11  13]  %%%%%%%%THE  INCOMPLETE  CASE  CSET%%%%%%%%%%%% 
%cset_rel=[4  5  6  7]%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o0/o0/o0/o0/o0/c%%°/o 

275  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%%%0/o0/o0/o0/o%0/o 

%%% 

%cset=[4  5  6  7  8  9  10  11]%%%%%%%%%%%%%THE  COMPLETE  CASE  CSET%%% 

%cset_rel=[4  5  6  7  8  9  10  1  1)%VoVoVo%%%Vo%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%0/o0/o%%%% 

280        %%% 

cset_size=length(cset); 
saveLL 
clear  L 

clear  kexta  kexto  mexta  mexto 
285        pack 

skyfull=zeros(nuindof,nuindof); 
counts  1  ; 

for  index  1 = 1  :numdof 
for  index2=indexl  :nunidof 
290  skyfull(  index  l,index2)=count; 

skyfull(index2,indexl  )=count; 
count=count+l; 
end 
end 
295        skyred=zeros(aset_size,aset_size); 
count=  1     ; 
for  indexl=l:aset_size 
for  index2= index  l:aset_size 

skyred(  index  1  ,index2  )=count; 
300  skyred(index2.  index  1  )=count, 
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count=count+l; 
end 
end 

skycset=zeros(cset_size,cset_size); 
305        count=l    ; 

for  indexl=l  :cset_size 
for  index2=indexl:cset_size 

skycset(  index  1  ,index2  )=counl; 
skycsel(index2,indexl  )=count; 
310  count=count+ 1 ; 

end 
end 

full_holder=zeros(  1  ,(numdof*(nuindof+ 1  ))/2); 
315        red_holder=zeros(  1  ,(aset_size*(aset_size+ 1  ))/2 ), 
cset_holder=zeros(  1  ,(cset_size*(cset_size+ 1  ))/2 ); 

if  meters 
\vaitbar_handle=waitbar(0,'Computing  Experimental  FRF'); 
320        else 

disp('Getting  experimental  FRF') 
end 

H_EXP=rj; 

for  count=  1  :numpoints 
325  if  meters 

waitbar(count/numpoints); 
end 

z_exp=k_exp+j*w(count)*c_exp-w(  count  )A2*m_exp; 
h_exp=inv(z_exp); 
330  h_exp=h_exp(aset,aset); 

skyindex=  1 ; 
for  index  1 = 1  :aset_size 

for  index2=mdexl:aset_size 
red_holder(skyindex)=h_exp(indexl,index2); 
335  skyindex=skyindex+l; 

end 
end 

H_EXP=[H_EXP  red_holder']; 
end 
340        if  meters 

close(  vvaitbarjiandle) 
end 

save  H_EXP  H_EXP 
clear  H_EXP 
345        pack 

if  meters 

\vaitbar_handle=waitbar(0,'Coniputing  Experimental  Impedence'); 
else 
350  disp('Getting  experimental  Impedence') 

end 

Z_EXP=[]; 

for  count=  1  :numpoints 
if  meters 
355  waitbar(count/numpoints); 

end 

z_exp=k_exp+j*\v(count)*c_exp-\v(count)A2*m_exp; 
h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
360  z_e\p=inv(h_exp); 
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skyindex=l  ; 

for  index  1  =  1  :aset_size 

for  iiidex2=indexl  :aset_size 
red_holder(skyindex)=z_exp(indexl,index2); 
365  skyindex=skyindex+l; 

end 
end 

Z_EXP=[Z_EXP  redjiolder"]; 
end 
370        if  meters 

close(\vaitbar_handle) 
end 

save  Z_EXP  Z_EXP 
clear  Z_EXP 
375        pack 

if  meters 

waitbar_handle=vvaitbar(0,'Computing  Reduced  Analytic  FRf' ); 
else 
380  dispCGettmg  Reduced  FRF') 

end 

HANAL_RED=0; 
for  count=  1  :numpoints 
if  meters 
385  \vaitbar(count/numpoints); 

end 

z_anal_red=kstat+j*w(coiuit)*cstat-w(count)A2*mstat; 
h_anal_red=inv(z_anal_red); 
skyindex=l  ; 
390  for  index  1 = 1  :aset_size 

for  index2=indexl  :aset_size 
red_holder(skyindex)=h_anal_red(indexl,index2); 
skyindex=skyindex+ 1 ; 
end 
395  end 

H_ANAL_REI>[H_ANAL_RED  redjiolder1]; 
end 

if  meters 
close(  \va  i  tbar_handle) 
400        end 

save  H_ANAL_R  H_ANAL_RED 

clear  Z_ANAL_RED  H_ANAL_RED  kexta  kexto  mexla  mexto 

pack 

405        %%  compute  and  save  L  matrix  diagonal  as  a  function  of  frequency%%%%%%%% 
if  meters 

waitbar_handle=\vaitbar(0,'Compuling  L  Diagonals'); 
else 
disp('Getting  L  Matrix') 
410        end 

L_DIAGS=[]; 
for  count=  1  :numpoints 
if  meters 

waitbar(count/numpoints); 
415  end 

z_anal=k_anal+j*w(count)*c_anal-w(coitnt)A2*m_anal; 
z_anal_red=kstat+j*w(count)*cstat-w(coiuit)A2*mslat; 
h_anal_red=inv(z_anal_red); 
z_exp=k_exp+j  *c_exp-\v(count)A2  *m_exp; 
420  h_exp=mviz_exp); 
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h_exp=h_exp(aset,aset); 

L=z_anal_red*(h_anal_red-h_exp)*z_anal_red; 

temp_l_diags=diag(L); 

L_DIAGS=[L_DIAGS  temp_l_diags(:)]; 
425        end 

if  meters 

close(waitbar_handle) 
end 

save  L_DIAGS  L_DIAGS 
430        clear  LDIAGS  z_anal  temp_l_diags  L 

clear  zexp  z_anal_red  u  v  hexp  h_anal_red 
pack 

if  meters 
435  waitbar_handle=waitbar(0,'Computing  Z  Analytic'); 

else 

disp('Getting  Z  anal') 
end 

Z_ANAL=[]; 

440        for  count=  1  .numpoints 

if  meters 

waitbar(count/numpoints); 
end 

z_anal=k_anal+j*w(count)*c_anal-w(count)A2*m_anal; 
445  skyindex=l    ; 

for  index  1 = 1  :numdof 

for  index2=indexl:numdof 
full_holder(skyindex)=z_anal(indexl,mdex2); 
skyindex=skyindex+ 1 ; 
450  end 

end 

Z_ANAL=[Z_ANAL  full_holder*]; 
end 

if  meters 
455  close(waitbar_handle) 

end 

save  Z_ANAL  Z_ANAL 
clear  Z_ANAL  z_anal 
pack 
460        if  meters 

\vaitbar_handle=waitbar(0,'Computing  Reduced  Z  Analytic'); 
else 

disp('Getting  Reduced  Z  anal') 
end 
465        Z_ANAL_RED=[]; 

for  count=  1  :numpoints 
if  meters 

waitbar(count/numpoints); 
end 
470  z_anal_red=kstat+j  *  w(count)*cstat-\v(count)A2  *mstat; 

skyindex=l     ; 
for  index  1 = 1  :aset_size 

for  index2=indexl:aset_size 
red_holder(skyindex)=z_anal_red(indexl,index2); 
sky  index=skyindex+ 1 ; 
end 
end 

Z_ANAL_RED=[Z_ANAL_RED  redjiolder']; 
end 
if  meters 
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close(waitbarhandle) 
end 

save  Z_ANAL_R  Z_ANAL_RED 
clear  Z_ANAL_RED  z_anal_red 
485        pack 

%%%computeDZ%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o%%%0/o0/o%0/o% 
if  meters 
waitbar_handle=waitbar(0,'Computing  DZ. . .'); 
490        else 

disp('Getting  DZ') 
end 
DZ=[]; 

for  count=  1  :numpoints 
495  if  meters 

waitbar(count/numpoints); 
end 

z_anal_red=kstat+j  *  w(count)*cstat-w(  count  )A2  *mstat; 
h_anal_red=inv(z_anal_red); 
500  z_exp=k_exp+j  *  w(count)*c_exp-\v(count)A2  *m_exp; 

h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
hacc=h_anal_red(cset_rel,cset_rel); 
hxcc=h_exp(cset_rel,cset_rel); 
505  dz=inv(inv(inv(hacc)*(liacc-hxcc)*inv(hacc)>hacc); 

skyindex=l     ; 
for  index  1 = 1  :cset_size 

for  index2=indexl:cset_size 
cset_holder(skyindex)=dz(  index  1  ,index2 ); 
510  skyindex- skyindex+1; 

end 
end 

DZ=[DZ  csetjiolder"]; 
end 
515        if  meters 

close(waitbarhandle) 
end 

save  DZ  DZ 

clear  hacc  lixcc  dz  hexp  z_anal_red  h_anal_red  z_exp 
520        pack 

%%%decompose  DZ  mto  DM,  DK,  &  rc%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if  meters 
waitbar_handle=\vaitbar(0,'DecomposingDZ...'); 
525        else 

disp('Decomposing  DZ') 
end 
DK=[]; 
DM=G; 
530        DC=[]; 

dz  1  =zeros(cset_size,cset_size); 
dz2=zeros(cset_size,cset_size); 
dz3=zeros(cset_size,cset_size); 
for  i=l  :numpoints-2 
535  if  meters 

waitbar(i/(numpoints-2)) 
end 

dztemp=DZ(ndx3d([l  cset_size*(cset_size+ 1  )/2  numpoints],l ,  1  :cset_size*(cset_size+l )/2,i)); 
skyindex=  1 ; 
540  for  mdexl  =  l:cset  size 
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for  index2=indexl:cset_size 
dz  1  (index  1  ,index2  )=dztemp(  skyindex ); 
dzl  (index2,index  1  )=dztemp(  sky  index); 
skyindex=skyindex+ 1 ; 
545  end 

end 

dztemp=DZ(ndx3d([  1  cset_size*(cset_size+ 1  )/2  munpoints],  1 , 1  :cset_size*(cset_size+ 1  )/2,i+ 1 )); 
skyindex=l; 
for  index  1 = 1  :cset_size 
550  forindex2=indexl:cset_size 

dz2(index  1  ,index2)=dztemp(  skyindex); 
dz2(index2,index  1  )=dztemp(skyindex); 
skyindex=skyindex+ 1 ; 
end 
555  end 

dztemp=DZ(ndx3d([  1  csel_size*(cset_size+ 1  )/2  numpoints],  1 , 1  :cset_size*(cset_size+ 1  )/2,i+2 )); 

skyindex=l; 

for  index  1 = 1  :cset_size 

for  index2=indexl  :cset_size 
560  dz3(indexl,index2)=dztemp(skyindex); 

dz3(  index2  ,index  1  )=dztemp(  skyindex); 
skyindex=skyindex+ 1 ; 
end 
end 
565  temp=[eye(cset_size)  -w(i)A2*eye(cset_size)   +j*w(i)*eye(cset_size); 

eye(cset_size)  -w(i+l  )A2*eye(cset_size)  +j*w(i+l  )*eye(cset_size); 
eye(cset_size) -w(i+2  )A2  *eye(cset_size) +j  *vv(  i+2  )*eye(cset_size)]\[dz  1  ;dz2;dz3]; 
ktemp=temp(  1  :cset_size,:); 
mtemp=temp(cset_size+l:2*cset_3ize,:); 
570  ctemp=temp(2*cset_size+l  :3*cset_size,:); 

skyindex- 1     ; 

for  index  1 = 1  :cset_size 

for  index2=indexl:cset_size 
575  cset_holder(  skyindex  )=ktenip(  index  1  ,index2); 

skyindex=skyindex+ 1 ; 
end 
end 

DK=[DK  cset_holder(:)]; 
580 

skyindex=l 

for  indexl = 1  :cset_size 

for  index2=indexl  :cset_size 
cset_holder(  skyindex  )=mtemp(  indexl, index2); 
585  skyindex=skyindex+l; 

end 
end 
DM=[DM  cset_holder(:)]; 

590  skyindex=l 

for  indexl = 1  :cset_size 

for  index2=indexl:cset_size 
cset_holder(skyindex)=ctemp(  indexl, index2); 
skyindex=skyindex+l ; 
595  end 

end 

DC=[DC  cset_holder(:)]; 
end 

if  meters 
close(\vaitbar_handle) 
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end 

save  DM  DM 

save  DK  DK 

save  DC  DC 
605        clear  DC  DK  DM 

save  exp  cexp  k_exp  m_exp 

clear  c_exp  k_exp  mexp 

save  stat  kstat  mstat  cstat 

clear  kstat  mstat  cstat 
610        save  anal  k_anal  m_anal  c_anal 

clear  k_anal  m_anal  c_anal 

pack 

if  complete 
615  loadZ_ANAL 

Z=Z_ANAL; 
clear  Z_ANAL 
else 
loadZ_ANAL_R 
620  Z=Z_ANAL_RED; 

clear  Z_ANAL_R£D 
end 

CORRHA=[]; 

tempza=zeros(aset_size,aset_size); 
625        tempz=zeros(cset_size,cset_size); 
if  meters 

waitbar_handle=waitbar(0,'lnstalling  DZ. . .'); 
else 
disp('Installing  DZ) 
630        end 

for  i=  1  :numpoints 
if  meters 

waitbar(  i/numpoints) 
end 
635  ztemp=Z(ndx3d([  1  aset_size*(aset_size+ 1  )/2  (numpoints-(i- 1 ))],  1 , 1  :aset_size*(aset_size+ 1  )/2, 1 )); 

skyindex=l; 
for  index  1  =  1  :aset_size 

for  index2=indexl:aset_size 
tempza(  index  1  ,index2  )=ztemp(  sky  index); 
640  tempza(index2, index  1  )=ztemp(skyindex); 

skyindex=skyindex+l ; 
end 
end 

zciMl; 

645  ztemp=DZ(ndx3d([  1  cset_size*(cset_size+ 1  )/2  (numpoints-<i- 1 ))],  1 , 1  :cset_size*(cset_size+ 1  )/2, 1 )); 

skyindex=l; 
for  index  1 = 1  :cset_size 

for  index2=indexl  xsetsize 
tempz(  index  1  ,index2)=ztemp(skyindex); 
650  tempz(mdex2,index  1  )=ztemp(skyindex); 

skyindex=skyindex+ 1 ; 
end 
end 

tempza(cset_rel,cset_rel)=tempza(cset_rel,cset_rel)+tempz; 
655  DZ(:,1)=[]; 

tempha=inv(tempza); 
CORRHA=[CORRHA  tempha(:)]; 
end 

if  meters 
660  close(waitbar_handle) 
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end         ' 

clear  DZ  Z 

save  CORRHA  CORRHA 

clear  CORRHA 

load  DK; 
load  DM; 
load  DC; 

670        load  slat 

CORRHAD=[]; 

tempza=zeros(aset_size,aset_size); 
tempz=zeros(cset_size,cset_size); 
tempk=zeros(cset_size,cset_size); 
675        tempm=zeros(cset_size,cset_size); 
tempc=zeros(cset_size,cset_size); 
if  meters 

waitbarjiandle=vvaitbar(0,'lnstalling  DK,  DM,  &  DC...'); 
else 
680  dispCIiistalliug  DK,  DM,  DC) 

end 

for  i=  1  :niunpoints-2 
if  meters 

waitbar(i/(numpoints-2)) 
685  end 

kcorrected=kstat; 

mcorrected=mstat; 

ccorrected=cstat; 

690  temp=DK(ndx3d([l  cset_size*(csel_size+ 1  )/2  (numpoints-(2))],l,l:cset_size*(cset_size+l)/2,l)); 

skyindex=l; 
for  index  1 = 1  xsetsize 

for  index2=indexl:cset_size 
tempk(  index  1  ,index2  )=temp(skyindex); 
695  tempk(index2, index  1  )=temp(  skyindex); 

skyindex=skyindex+l ; 
end 
end 

700  DK(:,1)=[]; 

temp=DM(ndx3d([  1  cset_size*(cset_size+ 1  )/2  (numpoints-(2))],  1 , 1  :cset_size*(cset_size+l  )/2, 1 )); 

skyindex=l; 

for  index  1  =  1  :cset_size 

for  index2=mdexl:cset_size 
tempm(indexl,index2)=temp(skyindex); 
tempm(  index2  .index  1  )=temp(  skyindex); 
skyindex=skyindex+ 1 ; 
end 
end 
710 

DM(:,1M]; 

temp=DC(ndx3d([  1  cset_size*(cset_size+l  )/2  (numpoints-(2))],  1 , 1  :cset_size*(cset_size+l  )/2, 1 )); 
skyindex=l; 
for  indexl=l:cset_size 
715  for  index2=index  1  :cset_size 

tempc(indexl,index2)=temp(  skyindex); 

tempc(index2,indexl)=teinp(  skyindex); 

skyindex=skyindex+l ; 
end 
end 
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DC(:,1)=[]; 

kcon-ected(cset_rel,cset_rel)=kcorrected(cset_rel,cset_rel)+tenipk; 
mcorrected(cset_rel,cset_rel)=incoiTected(cset_rel,cse(_rel)+teinpin; 
ccorrected(cset_rel,cset_rel)=ccorrected(cset_rel,cset_rel)+teinpc; 


tempza=kcorrected+j*w(i+l  )*ccorrected-vv(i+l  )A2*mcorrected; 

tempha=inv(tempza); 

CORRHAD=[CORRHAD  tempha(:)]; 
end 
730        if  meters 

close(vvaitbar_handle) 
end 

clear  tempk  tempm  tempc  tempza  tempha 
clear  kcorrected  mcorrected  ccorrecled 
735        save  CORRHAD  CORRHAD 
clear  CORRHAD  DK  DM  DC 

%%%save  the 

Workspace%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%0/o0/o0/o0/o%%%%0/o% 
740        save  INT 
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clear 
closeall 
load  INT 
titles=0; 
5        if  pswitch==y 

whitebgCwhite') 
close 
end 
h=0; 
10        if  complete 
fignum=  1 ; 
else 

fignum=21; 
end 
1 5        %%  write  frequencies  and  system  diag  to  diary  file  prt.out%%%%%%%%%%%%%%%%% 
ifexist('figOOO.ouf) 
delete  figOOO.out 
end 

diary  figOOO.out 
20        fprintfC     \ii') 
fprintfC     \n') 
fprintf('A  set\n') 
for  index=l:length(aset) 
fprintf([blanks(4-lengtli(int2str(aset(index)))),int2str(aset(  index))]) 
25  '         if  rem(index,12)==0 
fprintf(V) 
end 
end 

fprintf(V) 
30        fprintf(V) 

fprintf('0  set\n*) 
for  index=l:lengtli(oset) 
fprintf([blanks(4-lengm(int2str(oset(index)))),int2str(oset(  index))]) 
if  rem(  index,  12)==0 
35  fprintf(V) 

end 
end 

fprintf(*\n') 
fprintfCV) 
40        fprintf('Computed  C  sef\ri) 
for  index=  1  :length(cset) 
fpriiitf([blaiiks(4-lei'igtli(int2str(cset(index)))),int2str(cset(index))]) 
if  rem(  index,  1 2  )==0 
fprintf(V) 
45  end 

end 

fprintfCV) 
fprintfCV) 
fprintf(True  C  set\n') 
50        for  index=  1  :length(true_cset) 

fpriiitf([blanks(4-lengdi(int2str(true_cset(index)))),int2str(true_cset(mdex))]) 
if  rem(index,  12)=0 

fprintf(V) 
end 
55        end 

fprintf(V) 
fprintf(V) 

bigger=maxClength(omegaa),length(omegax)); 
if  complete 
60  allfreqs=[ ... 
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[omegaa/(2*pi);  zeros(bigger-length(omegaa),l)] ... 
(omegax/(2*pi);  zeros(bigger-length(omegax),  1 )] ... 

]; 

else 
65  allfreqs=[ ... 

[omegaa/(2*pi);  zeros(bigger-length(omegaa),l)      ]... 
[omegax/(2*pi);  zeros(bigger-length(omeg&x),l)        ]... 
[omegared/(2*pi);  zeros(bigger-length(omegared),  1 )  ]  ... 
[omegaexto/(2*pi);  zeros(bigger-length(  omegaexto),  1 )]... 
70  ]; 

end 

[aaa  bbb]=size(allfreqs); 
fprintf('System  Frequencies  (Hz)\n') 
if  complete 
75  fprintfC      Anal     Exp\ri) 

bbb=2; 
else 
fprintfC    Anal      Exp      Red     Oservn1) 
bbb=4; 
80        end 

for  indexl=l:aaa 
prtline=[]; 
forindex2=l:bbb 
if  allfreqs(  index  1  ,index2  )==0 
85  prtfreq-         '; 

else 
pit  freq=[blanks(8-lengtli(sprintf('%4g',allfreqs(indexl,index2 )))),... 
sprintfC%4g',allfreqs(  index  1  ,index2 ))] ; 
end 
90  prtline=[prtline,  prtfreqj; 

end 

prtline=[prtline,  WJ; 
fprintf(prtline) 
end 
95        diary  off 

%dos(type  figOOO.out  >  lptl  &'); 

xx=omegax(find(omegaa(lowmode)  <=  omegax  &  omegax  <=  omegaa(highmode)))/(2*pi); 
yx  1  =ones(length(xx),  1 ); 

aa=omegaa(find(omegaa(lowmode)  <=  omegaa  &  omegaa  <=  omegaa(highmode)))/(2*pi); 
100        yal=ones(lengm(aa),l); 

red=omegared(find(omegaa(lowmode)  <=  omegared  &  omegared  <=  omegaa(highmode)))/(2*pi); 
yred  1  =ones(  lengtli(red),  1 ); 

exto=omegaexto(find(omegaa(lowinode)  <=  omegaexto  &  omegaexto  <=  omegaa(lughniode)))/(2*pi); 
yexto  1  =ones(length(exto),  1 ); 
1 05        if  complete 

xO=max([xx(length(xx))aa(length(aa))]); 
else 

xO=max([xx(length(xx))  aa(lengtli(aa))  red(length(red))  exto(length(exto))]); 
end 
110        clear  omegax  omegaa  omegared  omegaexto 

clear  L  ZANAL  canal  cexp  conn  hexp  h_anal_red  kanal  kexp 
clear  kexta  kstat  kexto  m_anal  m_exp  mexto  mstat  temp_l_diags 
if  complete 
mid_index=round(length(cset)/2); 
115        else 

nud_index=roimd(lengtli(cset)/2); 
end 

load  H_ANAL_R 
loadH_EXP 
120 
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h=figure(h+l); 

fignuni=fignuin+l  ; 

plot(w/(2*pi),loglO(abs(H_ANAL_RED(nA\3d([l  aset_size*(asel_size+l  )/2  manpoinls], 

l,sk>Ted(cset_reI(mid_index),cset_rel(mid_index)),':')))),'r-',... 
1 25  w/(2*pi),loglO(abs(H_EXP(ndx3d([  1  aset_size*(aset_size+ 1  )/2  numpoinls],... 

l,skyred(cset_rel(nud_index),cset_rel(niid_index)),':')))),'b-.') 
v=axis; 

ylabel([,H(',iiit2str(cset(niid_mdex)),',,,mt2slr(cset(iiiid_ijidex));)  in/lbf (Log  1 0  of)']) 
xlabel('Oniega  (Hz)') 
130        if  titles 

title('Analytic  FRF  vs  Experimental  FRF') 
end 

hold  on 

ifabs(v(4)-v(3))<=  1 
135  v2(l)=v(l); 

v2(2)=v(2); 

v2(3)=v(3>.l*abs(v(4)); 

v2(4)=v(4>+.l*abs(v(4)); 

axis(v2); 
140        end 

grid  on 

lui=slegend('ANALYnC','EXPERIMENTAL'); 
axes(lili); 
hold  off 
145        ifpswitch=='y' 

prtfig(flgnuin) 

delete(h) 
end 

1 50        clear  H_ANAL_RED 
clear  H_EXP 

if  complete 
h=figure(h+l); 
155  fignum=fignum+ 1 ; 

else 
h=figure(h+2); 
fignum=fignum+2; 
end 
160        loadL 

plot(aset,abs(diag(L)),aset,abs(diag(L)),'*,) 
ylabel('  L(DOF)  lbf/in') 
if  complete 
xlabel('DOF') 
165        else 

xlabel('ASET  DOF) 
end 
grid  on 
if  titles 
1 70  title( ['Localization  Matrix  Diagonal  at  Omega  =  ',num2str(\vfreq/(2*pi)), '  Hz']); 

end 

if  pswitch==y 
prtfig(fignum) 
delete(h) 
175        end 

clear  L 

h=figure(h+l); 
fignum=fignum+ 1 ; 
180        loadL  DIAGS 
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mesh(w/(2*pi),aset,loglO(abs(L_DIAGS))) 
if  titles 

title('Frequeiicy  Dependence  of  Localization  Matrix  Diagonals') 
end 
185        xlabel('Omega  (Hz)') 
ylabelCDOF) 

zlabelCUDOF,  omega)  lbf/m  (Log  10  of)') 
grid  on 

ifpswitch— y 
1 90  prtfig(fignuin) 

delete(h) 
end 

h=figure(h+l); 
195        fignum=fignum+ 1 ; 

subplol(2,l,l) 

plot(w/(2*pi),loglO(L_DIAGS(ndx3d([aset_size  1  nunipoints],cset_rel(niid_index+2),l, ':')))) 

ylabel(  ["Error  Coord  L(',int2str(cset(imd_index+2)),',,,int2str(cset(inid_index+2 )),')']) 

if  titles 
200  title(*lbf/in  (Log  1 0  of)' ); 

end 

grid  on 

vl=axis; 

subplot(2,l,2) 


205 


210 


if  length(cset_rel)  <  length(aset) 
holdset=l:length(aset) ; 
holdset(cset_rel)=[] ; 
index=round(length(holdset)/2); 


plot(w/(2*pi),loglO(L_DIAGS(ndx3d([aset_size  1  numpoints],holdset(  index),  1, ':')))) 
ylabel(  ["Non-Error  Coord  L(',mt2str(aset(holdset(index))),',',int2str(aset(holdset(  index))),')']) 
if  titles 
titleCFrequency  Dependence  of  L  Matrix  Non-Error  set  Diagonal  Elements' ) 
215  end 

v=axis; 
hold  on 
v2(l)=v(l); 
v2(2)=v(2); 
220  v2(3)=v(3)+((v(4>v(3))*.5)-(vl(4>vl(3))*.5 ; 

v2(4)=v(3M(v(4)-v(3))*.5>f(vl(4)-vl(3))*.5; 
axis(v2), 
grid  on 
end 
225        xlabel('Omega  (Hz)') 
hold  off 
if  pswitch=='y' 
prtfig(fignum) 
delete(h) 
230        end 

clear  L_DIAGS 

load  Z_ANAL_R 

load  Z_EXP 
235        h=figure(h+l); 

fignum=fignum+ 1 ; 

plot(w/(2*pi),logl0(Z_ANAL_RED(ndx3d([l  aset_size*(asel_size+I  )/2 

iuunpoinls],l,skyred(cset_rel(nud_mdex),cset_rel(nud_index)),':'))),'r-',... 
W(2*pi),loglO(Z_EXP(ndx3d([l  aset_size*(aset_size+l  )/2 
240        numpoints],  1  ,sk>Ted(cset_reI(inid_index),cset_rel(mid_nidex)),':'))),'b-.') 
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xlabel('OmegaA2  (Hz)') 

ylabel([,Z(',mt2st^(cset(nlid_mdex)),',',int2str(cset(nlid_index)),,)  lbf/in  (Log  10  of)']) 
if  titles 

title('Aiialytic  Impedance  vs  Experimental  Impedance') 
245        end 

v=axis; 
hold  on 
ifabs(v(4>v(3))<=  1 

v2(l)=v(l); 
250  v2(2)=v(2); 

v2(3)=v(3>.l*abs(v(4)); 

v2(4)=v(4)Kl*abs(v(4)); 

axis(v2); 
end 
255        grid  on 

if  complete==  0 

hh=slegend('ANALYTIC  IMPEDANCE','EXPERIMENTAL  IMPEDANCE'); 
else 

hh=slegend('ANALYTIC  IMPEDANCE',EXPERIMENTAL  IMPEDANCE'); 
260        end 

axes(hh) 
hold  ofT 
ifpswitch==y 

prtfig(fignum) 
265  delete(h) 

end 
clear  Z_ANAL_RED  Z_EXP 

load  true_errs 
270        load  DK 

h=figure(h+l); 
fignum=fignum+l  ; 
if  complete 
plot(w(  1  :numpoints-2)/(2*pi),DK(ndx3d([  1  cset_size*(cset_size+ 1  )/2  numpoints-2],  1  ,skycset(mid  index,mid_index),':')),'r- 
275 

w(l:nunipoints-2)/(2*pi),true_stiffiiess(cset(mid_index),cset(nnd_index)).*ones(l,nuinpomts-2),'b-.') 
ylabel([DK(',int2str(cset(mid_index)),',',int2str(cset(imd_index));)  lbf/iri]) 
xlabel('Omega  (Hz)') 
else 
280  plot(w(l:numpomts-2)/(2*pi),loglO(abs(DK(ndx3d([l  cset_size*(cset_size+l)/2  numpoints- 

2],l,skycsel(nud_index,mid_index),':')))),'r-',... 

w(l:numpoints-2)/(2*pi),logl0((true_stiffness(cset(mid_index),cset(mid_index))==  0)+... 
true_stiffness(cset(mid_mdex),cset(nud_index))).*ones(l,numpoints-2),'b-.') 
ylabel(['DK(',mt2str(cset(imd_index)),',',int2str(cset(imd_index)),')  Ibf/in  (LoglO  of)']) 
285  xlabelCOmega  (Hz)') 

end 

v=axis; 
if  titles 
title('Computed  Stiffness  Error  vs  True  Stiffness  Error") 
290        end 

v2(l)=v(l); 
v2(2)=v(2); 
if  complete 
ifabs(v(4>v(3))<l 
295  v2(3)=v(3)-100*abs(v(4)-v(3)); 

v2(4)=v(4)+100*abs(v(4>v(3)); 
v(3)=v2(3); 
v(4)=v2(4); 
300  end 
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if  abs(tnie_stiffness(cset(mid_index),cset(mid_index))  -  v(3))  <  .25*abs(v(4)-v(3)) 
v2(3)=v(3)-.5*abs(v(4)-v(3)) ; 
v2(4)=v(4); 
elseif  abs(tme_stiffness(cset(nud_index),cset(niid_indcx))  -  v(4))  <  .25*abs(v(4)-v(3)) 
305  v2(4)=v(4)+.5*abs(v(4)-v(3)) ; 

v2(3)=v(3); 
else 
v2(3)=v(3); 
v2(4)=v(4); 
310  end 

else 
ifabs(loglO((true_stiffness(cset(niid_index),cset(nud_index))  —  0)+... 
true_stifmess(cset(inid_index),cset(inid_index)))  -  v(3))  <  .25*abs(v(4)-v(3)) 
v2(3)=v(3)-.5*abs(v(4)-v(3)) ; 
315  v2(4)=v(4); 

elseif  abs(loglO((tnie_stiffiiess(cset(nud_index),cset(mid_index))==  0)+... 
true_stifmess(cset(mid_index),cset(mid_index)))  -  v(4))  <  .25*abs(v(4)-v(3)) 
v2(4)=v(4)+.5*abs(v(4)-v(3)) ; 
v2(3)=v(3); 
320  else 

v2(3)=v(3); 
v2(4)=v(4), 
end 
end 
325        axis(v2) 
grid  on 

hh=slegend('COMPUTED  STIFFNESS  ERROR',TRUE  STIFFNESS  ERROR'); 
axes(hh) 
if  pswitch=='y' 
330  prtfig(fignum) 

delete(h) 
end 

clear  DK 
load  DM 
335       h=figure(h+l)  ; 

fignuin=fignuin+l  ; 
if  complete 
plot(\v(l:nunipoints-2)/(2*pi),DM(ndx3d([l  cset_size*(cset_size+l)/2  nuinpoints-2],l,skycset(niid_index,inid_index),':')),'r- 

340  \v(l:nmnpomts-2)/(2*pi),true_mass(cset(nud_index),cset(mid_index)).*ones(l,nuinpoints-2),'b-.') 

ylabel([DM(',inastr(cset(mid_index)),',',int2str(cset(mid_index)),')  lbf-sec^/m']) 
xlabel('Omega  (Hz)') 
else 
plot(w(  1  :nuinpoints-2)/(2*pi),logl0(abs(DM(ndx3d([  1  cset_size*(cset_size+l  )/2  numpoints- 
345        2],  1  ,skycsel(nud_index,niid_index),':')))),'r— ',... 
vv(  1  :iiuinpoints-2  )/(2  *pi  ),log  1 0( . . . 

true_mass(cset(mid_index),cset(mid_index))  ==  0  ... 
)+true_mass(cset(mid_index),cset(mid_index)). . . 
350  ).  *ones(  1  .numpoints-^XtK') 

ylabel([T>M(',int2str(cset(mid_index)),',',inaslr(cse((mid_index)),')  lbf-secA2/in  (LoglO  o0']) 
xlabel('Omega  (Hz)') 
end 
v=axis; 
355        if  titles 

title('Computed  Mass  Error  vs  True  Mass  Error') 
end 

hold  on 
v2(l)=v(l); 
360        v2(2)=v(2); 
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if  complete 
ifabs(v(4>v(3))<l 
v2(3)=v(3>10000*abs(v(4)-v(3)); 
v2(4)=v(4)+10000*abs(v(4>v(3)); 
365  v(3)=v2(3); 

v(4)=v2(4); 
end 

if  abs(true_mass(cset(mid_index),cset(inid_index))  -  v(3))  <  .25*abs(v(4)-v(3)) 
v2(3)=v(3>.5*abs(v(4)-v(3)); 
370  v2(4)=v(4); 

elseif  abs(true_mass(cset(mid_mdex),cset(niid_mdex))  -  v(4))  <  .25*abs(v(4)-v(3)) 
v2(4)=v(4>f.5*abs(v(4>v(3)); 
v2(3)=v(3); 
else 
375  v2(3)=v(3); 

v2(4)=v(4); 
end 
else 
if  abs(loglO((true_mass(cset(mid_index),cset(mid_index))  ==  0)+... 
380  tnie_iuass(cset(mid_index),cset(mid_index)))  -  v(3))  <  .25*abs(v(4)-v(3)) 

v2(3)=v(3>.5*abs(v(4>v(3)) ; 
v2(4)=v(4); 
elseif  abs(log  1 0((true_mass(cset(nud_index),cset(nud_index))  =  0)+. . . 
lrue_mass(cset(nud_index))cset(inid_index)))  -  v(4))  <  .25*abs(v(4)-v(3)) 
385  v2(4)=v(4>+.5*abs(v(4)-v(3)) ; 

v2(3)=v(3); 
else 
v2(4)=v(4); 
v2(3)=v(3); 
390  end 

end 

axis(v2) 
grid  on 

lm=slegend('COMPUTED  MASS  ERROR',TRUE  MASS  ERROR'); 
395        axes(hh) 

ifpswitch='y' 
prtfig(fignum) 
delete(h) 
end 
400        clear  DM 

load  DC 
h=figure(h+l); 
llgnum=fignum+ 1 ; 
405        if  complete 

plot(w(l:numpoints-2)/(2*pi),imag(DC(ndx3d([l  cset_size*(cset_size+l)/2  numpoints- 
2],  1  ,skycset(mid_index>mid_mdex),':'))),'^-,,. . . 

w(  1  :numpoints-2)/(2*pi),imag(true_damping(cset(mid_index),cset(nud_index)).*ones(  1  ,numpoints-2)),'b-.') 
ylabel([T>C(',int2str(cset(mid_index)),V,int2str(cset(mid_index)),')  lbf-sec/in']) 
410  xlabel('Oniega  (Hz)') 

else 

plot(w(l:numpoints-2)/(2*pi),loglO(abs(DC(ndx3d([l  cset_size*(csel_size+l)/2  numpoints- 
2],  1  ,skycset(mid_index,mid_index),':')))),'r-',. . . 
w(  1  :numpoints-2  )/(2  *pi  ),log  1 0(abs( . . . 
415  (... 

true_damping(cset(mid_index),cset(mid_index))  =  0  ... 
>i-true_daniping(cset(mid_index),cset(mid_index)). . . 
).  *ones(  1  ^iumpoints-2  )),"b-.') 

ylabel([T)C(',int2str(cset(mid_mdex)),V,int2str(cset(mid_index)),')  lbf-sec/in  (LoglO  of)']) 
420  xlabel('Omega  (Hz)') 
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end 
v=axis; 
if  titles 
title('Computed  Damping  Error  vs  True  Damping  Error') 
425        end 

v2(l)=v(l); 
v2(2)=v(2); 
if  complete 
ifabs(v(4>v(3))<l 
430  v2(3)=v(3)-100*abs(v(4>v(3)); 

v2(4)=v(4)+100*abs(v(4)-v(3)); 
v(3)=v2(3); 
v(4)=v2(4); 
end 
435  if  abs(lrue_damping(cset(mid_index),cset(mid_index))  -  v(3))  <  .25*abs(v(4)-v(3)) 

v2(3)=v(3)-.25*abs(v(4)-v(3)) ; 
v2(4)=v(4); 
elseif  abs(trae_dampmg(cset(inid_index),cset(nud_index))  -  v(4))  <  .25*abs(v(4)-v(3)) 
v2(4)=v(4>f  .25*abs(v(4>v(3)) ; 
440  v2(3)=v(3); 

else 
v2(3)=v(3); 
v2(4)=v(4); 
end 
445        else 

if  abs(logl  0((true_damping(cset(mid_index),cset(mid_index))  ==  0)+. . . 
true_damping(cset(mid_index),cset(nud_index)))  -  v(3))  <  .l*abs(v(4>v(3)) 
v2(3)=(v(3>.25*abs(v(4)-v(3))) ; 
v2(4)=v(4); 
450  elseif  abs(loglO((true_damping(cset(mid_index),cset(mid_index))  —  0)+... 

true_damping(cset(nud_index),cset(nud_index)))  -  v(4))  <  .l*abs(v(4)-v(3)) 
v2(4)=(v(4>+.25*abs(v(4>v(3))) ; 
v2(3)=v(3); 
else 
455  v2(3)=v(3); 

v2(4)=v(4); 
end 
end 

axis(v2); 
460        grid  on 

hh=slegend('COMPUTED  DAMPENG  ERROR',TRUE  DAMPING  ERROR'); 
axes(hh); 
if  ps\vitch==y 
prtfig(fignum) 
465  delete(h) 

end 
clear  DC  truedamping  truemass  lrue_stiffness 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
470        %%%% 

loadH_EXP 

load  CORRHA 

%for  nud_index=  1  :min(cset_size,multplot) 
475        h=figure<h+l); 

figinun=fignum+ 1 ; 

plot(w/(2*pi),loglO(abs(H_EXP(ndx3d([l  aset_size*(aset_size+l )/2 

numpouUs],Lskyred(cset_rel(nud_index),cset_rel(mid_index)),':')))),'r~',... 
w/(2*pi)JoglO(abs(CORRHA(udx3d([aset_sizeaset_sizenumpoiiUs],cset_rel(mid_index),cset_rel(mid_index),':')))),'b-.') 
480  ylabel([,H(',mt2suicset(nud_mdex)),',',iiit2str(cset(nud_index)),')  in/lbf  (LoglO  o0']) 
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xlabelOOmega  (Hz)') 
if  titles 

litleCExperimental  FRF  vs  DZ  Corrected  FRF') 
end 
485        grid  on 
v=axis; 

Wi=slegendCEXPERIMENTAL  FRF/CORRECTED  FRF); 
axes(hh)     ; 
if  pswiteh— y 
490  prtfig(fignuin) 

delete(h) 
end 

clear  CORRHA 
495        load  CORRHAD 

h=ligure(h+l); 

figiium=fignuin+ 1 ; 

plot(w/(2*pi),loglO(abs(H_EXP(ndx3d([laset_size*(aset_size+l)/2 

numpoints],  1  ,skyred(cset_rel(mid_index),cset_rel(mid_index)),':')))),,r— ',. . . 
500  w(l:nuinpomts-2)/(2*pi),logl0(abs(CORRHAD(ndx3d([aset_size  aset_size  (numpoints- 

2)],cset_rel(mid_iiidex),cset_rel(nud_mdex),,:')))),,b-.') 

ylabel([,H(,,int2str(cset(mid_index)),',',int2str(cset(mid_index)),,)  in/lbf  (LoglO  of)']) 

xlabel('Omega  (Hz)') 

if  titles 
505  titleOExperimental  FRF  vs  DK/D^M)C  Corrected  FRF') 

end 

grid  on 

v=axis; 

Wi=slegend(EXPERIMENTAL  FRF','CORRECTED  ANALYTIC  FRF'); 
510        axes(hh); 

ifpswitch='y' 
prtfig(fignuni) 
delete(h) 

end 
515        clear  CORRHAD 

save  plotnum  fignum 
clear 


103 


CHAP3.M 


%%%%SSTCONF  FILE  FOR  A  SPATIALLY  COMPLETE  BEAM 
beammdl  ;%beammdl  is  a  cantilevered  20  element  beam  model 

%simmdl  ;%simmdl  is  a  masses  and  springs  model 

pswitch- n'       ;%%do  we  print? 
5        titles=l  ;%%display  titles 

meters=  1  ;%use  progress  meters 

whitebgCblack')  ;%switch  to  black  figure  background 
close 

Iowmode=l; 
10        highmode=10; 

static=  1  ;%%%%%%%%%%%%%reduction  mediod  0=Guyan,  1  =IRS,2=Extraction 
%%define  A  set  &  O  set%%%%%%%%%%%%%%%%%%%%%%%%%Vo%%%%%VoVo°/0 
asel=l:numdof; 

1 5        %oset=2:2:numdof;  %%A  SPATIALLY  INCOMPLETE  BEAM%%%%%%%%%% 

%oset=sort(oset);%%%%%%%%%%%%%%%%%%%%%%%%0/o%%%%%%%%0/o%%%%% 

oset=[];   %%%%A  SPATIALLY  COMPLETE  BEAM%%%%%%%%%%%%%%% 
20        save  sstconf 
sst 
plotsst 
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%%%%SSTCONF  FILE  FOR  A  SPATIALLY  INCOMPLETE  BEAM 
beaiiundl  ;%beanundl  is  a  caatilevered  20  element  beam  model 

%simmdl  ;%simmdl  is  a  masses  and  springs  model 

pswitch- n'      ;%%do  we  print? 
5        titles=l  ;%%display  titles 

meters=  1  ;%use  progress  meters 

whitebgCblack')  ;%switch  to  black  figure  background 
close 

lowmode=  1 ; 
10        highmode=10; 

static=  1  ;%%%%%%%%%%%%%reduction  meUiod  0=Guyan,  1  =IRS,2=Extraction 
%%define  A  set  &  O  set%%%%%%%%%%%%%%%%%%%%%0/o0/0%%0/0%%%%%%%°/0 
aset=  1  :numdof; 
1 5        oset=2:2:numdof;  %%A  SPATIALLY  INCOMPLETE  BEAM%%%%%%%%%% 

oset=sort(oset);%%%%%%%%%%%%%%%%0/o%%%%%%%%%%%%%%%%%%%%%% 

save  sstconf 
sst 
20        plotsst 

beammdl  ;%beammdl  is  a  caiitilevered  20  element  beam  model 

%simmdl  ;%simmdl  is  a  masses  and  springs  model 

pswitch- n'      ;%%do  we  print? 
25        titles=l  ;%%display  titles 

meters=  1  ;%use  progress  meters 

wliitebgCblack')  ;%switch  to  black  figure  background 

close 

lowmode=l; 
30        liiglunode=10; 

%%define  A  set  &  O  set%%%%%%%%%%%%%%%%%%%0/0%%%%%%%%%%%%%% 
aset=l:numdof; 

oset=2:2:numdof;  %%A  SPATIALLY  INCOMPLETE  BEAM%%%%%%%%%% 
35        oset=sort(oset);  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

static=2; 
save  sstconf 
dofig4_3 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%%°/o%%%%% 

%%%%%%%%%%%% 

%%%  COMPARE  SINGLE  MODE  MATRIX  AND  INTEGRAL  %%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%  SOLUTIONS  UNDER  VARIOUS  CONDITIONS       %%%%%%%%%%%%%%%%%%%%%%%%%%%% 

5  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%%%%0/o%%% 

1M«%%%%%%%%%%% 

clear 

closeall 

load  INT 
1 0        whitebgCwhite') 

clear  L  Z_ANAL  c_anal  c_exp  corui  h_exp  h_anal_red  k_anal  k_exp 

clear  kexta  kexto  m_anal  mexp  mexta  mextO  templdiags 

clear  true_damping  truemass  true_stiflhess  z_anal  z_anal_red  z_exp 

clear  L  Z_ANAL  c_anal  c_exp  conn  h_exp  h_anal_red  k_anal 
1 5        clear  kexta  kexto  m_anal  mexta  mextO  temp_l_diags 

load  H_EXP 
load  exp 

cset=[5  7  9  11  13  15] 
20        cset_rel=[3  4  5  6  7  8] 
cset_size=length(cset) 
mid_index=round(length(cset)/2); 
fineness=200 
ODZl=[]; 
25        tempi  =[]; 

odfreq=omegax(  1 ); 
center_freq=odfreq/(2  *pi) 
\vl=odfireq; 
odlength=2*pi; 
30        odivisions=2; 

lowerfreq=odireq-l  5  *odlengtli 
upperfreq=odfreq+2  5  *odlength 

ow  1  =odireq-.5  *odlengtli:odleiigtli/odi  visions:odfreq+.  5  *odleugtli; 
partition=owl/(2*pi) 
35        df=(odlength/odivisions)/(2*pi) 
for  couiit=  1  :lengtli(ovvl ) 

z_anal_red=kstat-owl(count)A2*mstat; 

h_aiial_red=inv(z_aiial_red); 

z_exp=k_exp+j  *owl  (count  )*c_exp-ovv  1  (count)A2  *m_exp; 
40  h_exp=inv(z_exp); 

h_exp=h_exp(aset,aset); 

hacc=h_aiial_red(cset_rel,cset_rel); 

lixcc=h_exp(cset_rel,cset_rel); 

dz=inv(inv(inv(hacc)*(hacc-lixcc)*mv(hacc))-hacc); 
45  ODZl=[ODZl;dz]; 

templ=[templ;  eye(cset_size)  -owl(couiit)A2*eye(cset_size)  +j*owl(count)*eye(csel_size)]; 

end 

%Uielefldz=ODZl 
50        %t!iebigone=templ 

DKDMDC 1  =temp  1  \ODZ  1 ; 

ODK 1  =DKDMDC  1  ( 1  :cset_size,: ) 

ODM 1  =DKDMDC  1  (cset_size+ 1 :2*cset_size,: ) 

ODC 1  =DKDMDC  1  (2*cset_size+ 1 :3*cset_size,:) 
55        save  ODM  1  ODM  1  DKDMDC  1  ODZ1  owl  wl 

saveODKl  ODK1 

save  ODC  1  ODC1 

%l  I  M  I  I  I  I  I  I  I  1  I  II  I  1  I  I  I  I  I  I  1  I  1  I  I  I   I   I   I  I   I   I   I  I  I  I  I  I   I   I  I  I   I  |  |   1   I   1   |   |   !   |  |  I  I  I   I  1  1   I  I  I  I   I  I  I  I  I  I 

60         load  stat; 
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odcorrha«=[]; 

w=lowerfreq:(upperfreq-lowerfreq)/fineness:upperfreq; 
numpoints=length(w); 
kcorrected=kstat; 
65        mcorrected=mstat; 
ccorrected=cstat; 

kcorrected(cset_rel,cset_rel)=kcorrected(cset_rel,cset_rel)+ODK  1 ; 
mcorrected(cset_rel,cset_rel)=mcorrected(cset_rel,cset_rel)+ODM  1 ; 
ccorrected(cset_rel,cset_rel)=ccorTected(cset_rel,cset_rel)+ODC  1 ; 
70        for  i=l:numpoints 

tempza=kcorrected+j*w(i)*cx;orrected-w(i)A2*mcorrected; 
tempha=inv(tempza); 
odcorrha=[odcorrha  tempha(:)]; 
end 
75 

uncorrha=[]; 
for  i=l:numpoints 
tempza=kstat+j*w(i)*cstat-w(i)A2*mstat; 
tempha=inv(tempza); 
80  uncorrha=[uncorrha  tempha(:)]; 

end 

if  meters 

waitbar_handle=waitbar(0,'Cornputing  Experimental  FRF'); 
85        else 

dispOGetting  experimental  FRF1) 
end 

H_EXP=[]; 
for  coimt=  1  :length(w) 
90  if  meters 

waitbar(coimt/length(w)); 
end 

z_exp=k_exp+j*w(count)*c_exp-w(count)A2*m_exp; 
h_exp=inv(z_exp); 
95  h_exp=h_exp(aset,aset); 

skyiiidex=  1 ; 
for  index  1 = 1  :aset_size 
for  index2=indexl:aset_size 
red_holder(skyindex)=h_exp(indexl,mdex2); 
100  skyindex=skyiiidex+l; 

end 
end 

H_EXP=[H_EXP  redjiolder*]; 
end 
105        if  meters 

close(waitbar_handle) 
end 
clear  tempha  tempza  kcorrected  mcorrected  ccorrected 

110        figure(l); 

plot(w/(2*pi),loglO(abs(H_EXP(ndx3d([l  aset_size*(aset_size+l )/2  numpoints],... 

1  ,skyred(cset_rel(mid_index),cset_rel(imd_iiidex)),':')))),'r-',. . . 
w/(2*pi)Jogl0(abs(odcorrha(ncbc3d([aset_sizeaset_sizenmnpomts],cset_rel(imd_mdex),cset_rel(nud_mdex),':,)))),'g-', 

w/(2*pi)Jogl0(abs(uria)nlia(ndx3d([aset_sizeaset_sizenumrx>mts],cset^ 

ylabel(['H(',uit2str(cset(mid_index)),'),,mt2str(cset(mid_mdex));)  in/lbf  (log  1 0  of)']) 

120        xlabel('Omega  (Hz)') 
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%titleCSingle  Mode  Corrected  FRF  (Mode  1 ) ') 

grid 

hh=slegendCExperimental  FRF,  "Mode  1  Corrected  FRF',  'Uncorrected  Analytic  FRF'); 

125        axes(lih) 

print  -dinfile  fig5_l 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%%%%% 
%%%%®e©i®SasaS®Z  into  DM,  DK,  &  DC    %%%%%%%%%0/o0/o0/o0/o%%%0/o0/o%%%%%0/o%%%%%%0/o0/o0/o%% 
1 30        %%%  USING  MATRIX  TECHNIQUE  %%%%%%%0/o%0/o%%%%%0/o%%%%0/o0/o%%%0/o%0/o0/o%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%0/o%%0/o%%% 

templ=[]; 
odlengtii=2*pi; 
135        owl=omegax(l)-5*odlengUi:(omegax(2)-omega\(l)+10*odlengtli)/14:omegax(2)+5*odlengUi; 

for  count=  1  :lengUi(owl ) 

z_anal_red=kstat-ow  1  (count  )A2  "mslat; 

h_anal_red=inv(z_anal_red); 
1 40  z_exp=k_exp+j  *owl  (count)*c_exp-owl  (count  )A2*m_exp; 

h_exp=inv(z_exp); 

h_exp=h_exp(aset,aset); 

hacc=h_anal_red(cset_rel,cset_rel); 

hxcc=h_exp(cset_rel,cset_rel); 
1 45  dz=inv(inv(inv(hacc)*(liacc-hxcc)*inv(hacc))-hacc); 

ODZl=[ODZl;dz]; 

tempi  =[templ;  eye(cset_size)  -owl  (count  )A2*eye(cset_size)  +j*o\vl(count)*eye(cset_size)]; 
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end 


DKDMDCl=templ\ODZl; 
ODK1  =DKDMDC  1  ( 1  :cset_size,:) ; 
ODMl=DKDMDCl(cset_size+l :2*cset_size,:) ; 
ODCl=DKDMDCl(2*cset_size+l:3*cset_size,:); 
155        save  ODM1  ODM1  DKDMDC1  ODZ1  owl  wl 
saveODKlODKl 
saveODCl  ODC1 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%°/o%%% 
160        UMm^/o%%%%%%% 
odcorrlia=[]; 

w=omegax(l)-15*pi:(oinegax(2)-omegax(l)+35*pi)/fineness:oinegax(2)+20*pi; 
numpoints=lengUi(w); 
kcorrected=kstat; 
165        mcorrected=mstat; 
ccorrected=cstat; 

kcorrected(cset_rel,cset_rel  )=kcorrected(cset_rel,cset_rel  )+ODK  1 ; 
mcorrected(cset_reLcset_rel  )=mcorrected(cset_rel,cset_rel  >K)DM  1 ; 
ccorrected(cset_rel,cset_rel)=ccorrected(cset_rel,cset_rel)K)DC  1 ; 
1 70        for  i=  1  :numpoints 

tempza=kcorrected+j*w(i)*ccorrected-w(i)A2*mcorrected; 
tempha=inv(tempza); 
odcorrha=[odcoiTha  tempha(:)]; 
end 
175 

uncorrha=[]; 
fori=l:numpoints 
tempza^kstat+j  *  w(  i)*  cstat-w(i  )A2  *mstat; 
tempha=inv(tempza); 
180  uncorrha=[uncorrha  temphaQ], 
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end 

if  meters 

waitbar_handle=waitbar(0,'Computing  Experimental  FRF'); 
185        else 

disp('Getting  experimental  FRF) 
end 

H_EXP=[]; 
for  count=  1  :length(w) 
190  if  meters 

waitbar(count/lengtli(w)); 
end 

z_exp=k_exp+j*\v(count)*c_exp-w(  count  )A2*m_exp; 
h_exp=inv(z_exp); 
1 95  h_exp=h_exp(aset,aset); 

skyindex=l; 
for  index  1 = 1  :aset_size 
for  index2=indexl:aset_size 
red_holder(skyindex)=h_exp(indcx  1  ,index2); 
200  skyindex=skyindex+l; 

end 
end 

H_EXP=[H_EXP  redjiolder1]; 
end 
205        if  meters 

close(\vaitbar_handle) 
end 
clear  tempha  tempza  kcorrected  mcorrected  ccorrected 

210        figure(2); 

plot(\v/(2*pi),logl0(abs(H_EXP(ndx3d([l  aset_size*(aset_size+ 1  )/2  numpomts],... 

1  ,skyred(cset_rel(mid_index),cset_rel(mid_index)),':')))),'r-',.  ■  • 

\v/(2*pi),logl0(abs(odcorrha(ndx3d([aset_size  aset_sizenumpoints])cset_rel(mid_index),cset_rel(mid_index),':')))),'g-',... 

\v/(2*pi),loglO(abs(uncorTha(ndx3d([aset_sizeaset_sizenumpoints],cset_rel(mid_index),cset_rel(mid_index),,:')))),'b:') 

ylabel(['H(',int2str(cset(mid_index)),',',int2str(cset(mid_index)),')  in/lbf  (loglO  of)']) 

220        xlabeK'Omega  (Hz)') 

%titleCSingle  Mode  Corrected  FRF  (Mode  1 ) ') 

grid 

lm=slegend(Experimental  FRF',  'Single  Mode  Corrected  FRF',  'Uncorrected  Analytic  FRF'); 

225        axes(hh) 

print  -dmfile  fig5_2 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%% 
%%%%«4&38f&d%mSZ  into  ODM,  ODK,  &  ODC  %%%%0/o%%%%0/o%%%%%%%%%0/o0/o%%%%%%%%0/o 
230        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%0/o0/o%%%% 
^tt^Mi%if^^(pp^f^eq-lo\verfreq)/fineness:upperfreq; 
lenctrl=[25  10  1]; 
odivisions=3; 
for  ind=l:length(lenctrl) 
235  odlengtli=lenctrl(ind)*2*pi; 

ODZ=[]; 
temp=[]; 

o%\'=odfreq-.5*odlengtli:odlengtli/odivisions:odfreq+.5*odlenglh; 
for  count=  1  :lengtl;(ow) 
240  z_anal_red=kstat-ow(  count  )A2*mslat; 
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h_anal_red=inv(z_anal_red); 

z_exp=k_exp+j  *ow(count)*c_exp-ovv(count)A2  *m_exp; 
h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
245  hacc=h_anal_red(cset_rel,cset_rel); 

hxcc=h_exp(cset_rel,cset_rel); 
dz=inv(inv(inv(hacc)*(hacc-lixcc)*inv(hacc))-hacc); 
ODZ=[ODZ;  dz]; 
temp=[temp;  eye(cset_size)  -ow(count)A2*eye(cset_size)  +j*ow(count)*eye(cset_size)]; 

end 

DKDMDC=temp\ODZ; 

ODK=DKDMDC(  1  :cset_size,:) ; 
255  ODM=DKDMDC(cset_size+l:2*cset_size,:) ; 

ODC=DKDMDC(2*cset_size+l:3*cset_size,:); 

clear  hacc  lixcc  h_exp  z_exp  z_anal_red  h_aiial_red  dz 

save  ODM  ODM  DKDMDC  ODZ  ow 

save  ODK  ODK 
260  save  ODC  ODC 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o0/o0/o0/o0/o%%0/o% 
%WftM6B4%%%% 
waitbar_handle=waitbar(0,,Computmg  Experimental  FRF); 
else 
265  disp('Getting  experimental  FRF) 

end 

H_EXP=[]; 

for  count=l  :lengtli(plotw) 
if  meters 
270  waitbar(count/length(plotw)); 

end 

z_exp=k_exp+j*plotw(count)*c_exp-plotw(count)A2*ni_exp; 
h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
275  skyindex=l; 

for  index  1  ■  1  :aset_size 
for  index2=indexl:aset_size 
red_holder(skyindex)=h_exp(  index  1  ,index2); 
sky  index=skyindex+ 1 ; 
280  end 

end 

H_EXP=[H_EXP  redjiolder']; 
end 

if  meters 
285  close(waitbar_handle) 

end 

odcorrha=[], 
kcorrected=kstat; 
290  mcorrected=mstat; 

ccorrected^cstat; 

kcorrected(cset_rel,cset_rel)=kcorrected(cset_rel,cset_rel)+ODK; 
mcorrected(cset_rel,cset_rel)=mcorrected(cset_rel,cset_rel)+ODM; 
ccorrected(cset_rel,cset_rel)=ccorrected(cset_rel,cset_rel)+ODC; 
295  for  i=  1  :lengtli(plot\v) 

tempza=kcon"ected+j*plotvv(i)*ccorTected-plotw(i)A2*mcorTected; 
tempha=inv(tempza); 
odcorrha=[odcorrha  tempha(:)]; 
end 
300  ifmd==l 
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plotl=odcorrha; 
elseif  iiid==2 

plot2=odcorrha; 
elseif  ind=3 
305  plot3=odcorrha; 

else 

plot4=odcorrha; 
end 
end 


310 


figure(3); 

plot(plolw/(2*pi),loglO(abs(H_EXP(ndx3d([  l  aset_size*(aset_size+l  )/2  length(plotw)],. 


l,skyred(cset_rel(mid_index),cset_rel(mid_index)),':')))),,r-',.-. 
31 5  plotw/(2*pi),  Iogl0(abs(plotl(ndx3d([aset_size  asel_size  lengui(plotw)],cset_rel(mid_index),cset_rel(mid_index),':')))),'g- 

plotw/(2*pi),  Iogl0(abs(plot2(ndx3d([aset_size  aset_size 
lengtli(plotw)],cset_rel(nud_index),csel_rel(inid_index),':')))),'b:',... 

plorw/(2*pi),  Iogl0(abs(plot3(ndx3d([aset_sizeaset_size  lengtli(plotw)],cset_rel(inid_index),cset_rel(mid_index),':')))),'k-.') 

ylabel(['H(',int2str(cset(imd_index)),V,int2str(cset(niid_index)),')  in/lbf(loglO  of)']) 

xlabelC'Omega  (Hz)') 
grid  on 
325        %UUe('EXPERMENTAL  FRF  AND  CORRECTED  FRFS  USING  MODE  5  SOLUTIONS') 
hh=legendCExp  FRF, ... 

[num2str(lenctrl(l))/  Hz  Bandwidth  Corrected  FRF'], ... 
[num2str(lenctrl(2)),'  Hz  Bandwiddi  Corrected  FRF*], ... 
[nuin2str(lenctrl(3)),'  Hz  BandwidUi  Corrected  FRF]) ; 
330        axes(hli) 

print  -dnifile  fig5_3 

odlengtli=l*2*pi; 
countctrl=[3  10  50  200] ; 
335        for  ind=  1  :leugtli(countctrl) 
odi  visions=countctrl(ind)- 1 ; 
ODZ=[]; 
temp=[]; 

ow=odfreq-.  5  *odlengUi:odlength/odi  visions:odfreq+.  5  *odlengtli; 
340  for  count=  1  :lengdi(ow) 

z_anal_red=kstat-ow(  count  )A2*ms  tat; 
h_anal_red=inv(z_anal_red); 

z_exp=k_exp+j  *ow(count)*c_exp-ow(  count  )A2  *m_exp; 
h_exp=inv(z_exp); 
345  h_exp=h_exp(aset,aset); 

hacc=h_anal_red(cset_rel,cset_rel); 
lixcc=h_exp(cset_rel,cset_rel ); 
dz=inv(inv(inv(hacc)*(hacc-lvxcc)*inv(hacc))-hacc); 
ODZ=[ODZ;  dz], 
350  temp=[temp;  eye(cset_size)  -ow(count)A2*eye(cset_size)  +j*o\v(counl)*eye(cset_size)]; 

end 

DKDMDC=temp\ODZ; 
355  ODK=DKDMDC(  1  :cset_size,:) ; 

ODM=DKDMDC(cset_size+l:2*cset_size,:) ; 

ODC=DKDMDC(2  *cset_size+ 1 :3  *cset_size,:); 

clear  hacc  hxcc  h_exp  z_exp  z_anal_red  h_anal_red  dz 

save  ODM  ODM  DKDMDC  ODZ  ow 
360  save  ODK  ODK 
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save  ODC  ODC 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% 
odcorrha=[]; 
365  kcorrected=kstat; 

mcorrected=mstat; 
ccorrected=cstat; 

kcorrected(cset_rel,cset_rel)=kcon'ected(cset_rel,cset_rel)+ODK; 
mcorrected(cset_rel,cset_rel)=nicx)iTected(cset_rel,cset_rel)+ODM: 
370  ccorrected(cset_rel,cset_rel)=ccorrected(cset_rel,cset_rel)+ODC; 

for  i=  1  :lengtli(plotw) 
tenipza=kcorrected+j  *plolw(i)*ccorrecled-plotw(i  )A2  "mcorrecled; 
lempha=inv(  tempza); 
odcorrha=[odcorrha  tempha(:)]; 
375  end 

if  ind==l 
plotl=odcorrha; 
elseif  ind—2 
plot2=odcorrha; 
380  elseif  ind=3 

plot3=odcorrha; 
else 

plot4=odcorrha; 
end 
385        end 

figure(4); 

plot(plot\v/(2*pi),loglO(abs(H_EXP(ndx3d([l  asel_size*(aset_size+l )/2  length(plotw)]:... 

390  l,skyred(cset_rel(mid_index),cset_rel(nud_index)),':')))),,r-',... 

plot\v/(2*pi),  Iogl0(abs(plotl(ndx3d([aset_size  aset_size  lengUi(plotw)],csetjel(inid_index),cset_rel(mid_index),':')))),'g- 

plotw/(2*pi),  Iogl0(abs(plot2(ndx3d([aset_size  aset_size 

lengtli(plot\v)],cset_rel(nud_index),cset_rel(mid_uidex),':')))),'b:',-- 
395  plotw/(2*pi),  logl  0(abs(plot3(ndx3d([aset_size  aset_size  length(plotw)],cset_rel(niid_iiidex),cset_rel(mid_index),':')))),'k- 

plot\v/(2*pi),  Iogl0(abs(plot4(ndx3d([aset_size  aset_size  length(plol\v)],cset_rel(inid_index),cset_rel(mid_index),':')))),'nv.') 

ylabel(['H(',int2str(cset(nud_index)),V,int2str(csel(nud_index));)  in/lbf  (loglO  o0']) 

xlabel('Omega  (Hz)') 
grid  on 

%titleCEXPERMENTAL  FRF  AND  CORRECTED  FRFS  USING  MODE  1  SOLUTIONS') 
hh=legend('Exp  FRF', ... 
405  [nuin2str(countctrl(  1 )),'  Pts  Corrected  FRF'], . . . 

[nuin2slr(countctrl(2)),'  Pts  Corrected  FRF], ... 

[nuni2str(countctrl(3)),'  Pts  Corrected  FRF'], ... 

[num2str(countctrl(4)),'  Pts  Corrected  FRF]); 
axes(hh) 
410        print  -dmfile  fig5_4 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%^as«%ffiwasm)z  into  dm,  dk,  &  dc  %%%%%%%%%%%%%o/o%%%%%%o/o%%%o/o%o/oo/o%%%o/o%%%% 

%%%  USING  INTEGRAL  TECHNIQUE    %%%%%%%%%%%%%%%0/o%%%0/o%0/o0/o0/o0/o%%%%0/o%%%0/o%0/o%%% 
415        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%plot\v=lowerfreq:dw:upperfreq; 
lenctrl=[25  10  1  ]; 
odivisions=3; 
420        forind=l:length(lenctrl) 
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odlength=lenctrl(ind)*2*pi; 

owl=odfreq-.5*odlengtli:odlengtli/(odivisions):odfreq+.5*odlengtli; 
Wl=owl; 

%%compute  integrals    %%%%%%%%0/o0/o0/o%%0/o%%%0/o%%0/o0/o%0/o%%0/o%%%0/o0/o%%0/o% 
425  %weiglit=ones(size(owl  ))*(omegax(  1 )); 

weight=owl; 
DZ_R=[]; 
DZ_I=[]; 

for  count=l  :lengtli(owl ) 
430  z_anal_red=kstat/(j*owl(count)>+-j*owl(count)*mstat; 

h_anal_red=inv(z_anal_red); 

z_exp=k_exp/(j  *owl(count))+c_exp/(j  *owl  (count))+j  *o\vl  (couiit)*m_exp; 

h_exp=inv(z_exp); 

h_exp=h_exp(aset,aset); 
435  hacc=h_anal_red(cset_rel,cset_rel); 

hxcc=h_exp(cset_rel,cset_rel); 

dz=iiiv(mv(uiv(liacc)*(hacc-hxcc)*inv(hacc))-hacc); 

dz_r=real(dz); 

dz_i=imag(dz); 
440  DZ_R=PZ_R  dz_r(:)]; 

DZ_I=[DZ_I  dz_i(:)]; 
end 

[INTK INTM  INTC]=intsub(DZ_I,DZ_R,owl  ,W  1 , 1  ./(j  *owl  ),cset_size); 
445  odcorrha=[]; 

load  stat 

kcorrected=kstat; 
mcorTected=nistat; 
ccorrected=cstat; 
450  kcorrected(cset_rel,cset_rel)=kcon"ected(cset_rel,cset_rel>i-INTK; 

mcoirected(cset_rel,cset_rel)=mcorrected(cset_rel,cset_rel)+INTM; 
ccorrected(cset_rel,cset_rel)=ccorrected(cset_rel,cset_rel)+INTC; 
for  i=l:length(plotw) 
tempza=kcorrected+j*plotw(i)*ccorrected-plotw(i)A2*mcorrected; 
455  tempha=inv(tempza); 

odcorrha=[odcorrha  tempha(:)]; 
end 

ifind==l 
intplot  1  =odcorrha; 
460  save  INTM1  INTM 

save  INTK1  INTK 
savelNTCl  INTC 
elseif  ind==2 
intplot2=odcorrha; 
465  save  INTM2  INTM 

save  INTO  INTK 
save  INTC2  INTC 
elseif  ind=3 
intploG=odcorrha; 
save  INTM3  INTM 
save  INTK3  INTK 
save  INTC3  INTC 
else 
intplot4=odcorrha; 
475  save  INTM4  INTM 

save  INTK4  INTK 
save  INTC4  INTC 
end 
end 
480 
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if  meters  ' 

waitbar_handle=waitbar(0,'Computing  Experimental  FRF); 
else 
disp('Getting  experimental  FRF') 
485        end 

H_EXP=[]; 

for  count=l  :length(plotw) 
if  meters 
waitbar(count/length(plotw)); 
490  end 

z_exp=k_exp+j  *plotw(count)*c_exp-plotw(count)A2  *m_exp; 
h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
skyindex=  1 ; 
495  for  index  1 = 1  :aset_size 

for  index2=indexl  :aset_size 
red_holder(skyindex)=h_exp(  index  1  ,index2); 
skyindex=skyindex+ 1 ; 
end 
500  end 

H_EXP=[H_EXP  redjiolder1]; 
end 

if  meters 
close(waitbar_handle) 
505        end 

figure(5); 

plot(plotw/(2*pi),loglO(abs(H_EXP(ndx3d([l  aset_size*(aset_size+l )/2  length(plotw)],... 

510  1  ,skyred(cset_rel(nud_index),cset_rel(mid_index)),,:')))),'r-',. . . 

plotw/(2*pi),  Iogl0(abs(intplotl(ndx3d([aset_sizeaset_size 
lengtli(plotw)],cset_rel(mid_index),cset_rel(nud_index),':')))),'g— ',... 

plotw/(2*pi),  Iogl0(abs(intplot2(ndx3d([aset_size  aset_size 
lengtli(plotw)],cset_rel(nud_index),cset_rel(mid_index),':')))),'b:',..- 
515  plotw/(2*pi),  Iogl0(abs(intplot3(ndx3d([aset_sizeaset_size 

length(plotw)],cset_rel(nud_index),cset_rel(nud_index),':')))),'k-.') 
ylabel([,H(',int2st^(cset(nud_index)),';,int2st^(cset(mid_index)),,)  in/lbf  (log  1 0  of)']) 

xlabel('Omega  (Hz)*) 
520        grid  on 

%titleCEXPERMENTAL  FRF  AND  CORRECTED  FRFS  USING  MODE  1  INT  SOLUTIONS') 
hh=legendOExp  FRF, ... 

[num2str(lenctrl(l)),'  Hz  Bandwidth  Corrected  FRF'],... 
[num2str(lenctrl(2)),'  Hz  Bandwidth  Corrected  FRF],... 
525  [num2str(lenctrl(3)),'  Hz  Bandwidth  Corrected  FRF']); 

axes(hh) 
print  -dmfile  fig5_5 

odlengtli=2*pi; 
530        countctrl=[3  10  50  200]; 
for  ind=  1  :lengtli(countctrl) 
odi  visions=countctrl(  ind); 

owl  =odfreq-.  5  *odlengtli:odlengtli/(odi  visions- 1  ):odfreq+.5*odlengtli; 
Wl=owl; 
535  %%compute  integrals     %%%%%%%%%%%%%%%%%%%%%0/o%%0/o%%%0/o%%%0/o%0/o%0/o%%%0/o%0/o0/o 

%weiglit=ones(size(owl  ))*(omegax(  1 )); 
weight=owl ; 
ifind==l 
owl 
540  df=(odlength/(udivisions-l  ))/(2*pi) 
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end 

%%  compute  integrals     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%o/0o/0%%%%%%%o/0o/0% 
DZ_R=[J; 
DZ_I=G; 
545  for  count=  1  :length(owl ) 

z_anal_red=kstat/(j  *owl  (count))+j  *owl  (couiit)*mslat; 
h_anal_red=inv(z_anal_red); 

z_exp=k_exp/(j  *owl  (count)>i-c_exp/(j  *owl  (count))+j  *owl  (count  )*m_exp; 
h_exp=inv(z_exp); 
550  h_exp=k_exp(aset,aset); 

hacc=h_anal_red(cset_rel,cset_rel); 
hxcc=h_exp(cset_rel,cset_rel); 
dz=inv(inv(inv(hacc)*(hacc-lixcc)*inv(hacc))-hacc); 
dz_r=real(dz); 
555  dz_i=imag(dz); 

DZ_R=[DZ_R  dz_r(:)]; 
DZ_I=[DZ_I  clz_i(:)]; 
end 

ifind==  1 
560  DZ_R; 

DZ_I; 

for  i=  1  :length(ow  1 ) 
dztempl=DZ_R(ndx3d([l  cset_size*(cset_size+l  )/2  length(owl  )],1 ,1  :cset_size*(cset_size+l  )/2,i)); 

565  dztemp2=DZ_I(ndx3d([  1  cset_size*(cset_size+ 1  )/2  lengtli(owl )],  1 , 1  :cset_size*(cset_size+ 1  )/2,i)); 

skyindex=l; 
for  index  1 = 1  :cset_size 
for  index2=indexl:cset_size 
570  dzl(indexl,index2):=dztempl(skyindex); 

dzl  (index2,index  1  )=dztemp  1  (skyindex); 
dz2(indexl  ,index2)=dztemp2(skyindex); 
dz2(index2,indexl)=dztemp2(skyindex); 
skyindex=skyindex+ 1 ; 
575  end 

end 
dzl 
dz2 
end 
580 

end 

[INTK INTM  INTC]=intsub(DZ_I,DZ_R,o\vl,Wl,l./(j*owl),cset_size); 
ifind==  1 

[MIX  INTM  INTC]=indsub(DZ_LDZ_R,owl ,  W 1 , 1  ./(j  *ow  1  ),cset_size); 
585  INTM 

INTK 
HNTC 
end 

590        %8lfla»«fi^/o%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o% 

°46tiabstat 

kcorrected=kstat; 

mcorrected=mstat; 

ccorrected^cstat; 
595  kcorTected(cset_rel,cset_rel)=kcorrected(cset_rel,cset_rel)+INTK; 

mcorrectedCcse^rekcse^rel^mcorrectedCcset^ekcse^reiyt-INTM; 

ccon"ected(cset_rel,cset_rel)=ccorrected(cset_rel,cset_rel)+INTC; 

for  i=l:lengtli(plotw) 
tempza=kcorrected+j*plotw(i)*ccorrected-plotw(i)A2*mcorTected; 
600  tempha=inv(tempza), 
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odcorrha=[odcorrha  tempha(:)]; 
end 

ifind==l 
iiitplot  1  =odcorrha; 
605  save  INTM1  INTM 

save  INTK1  INTK 
savelNTCl  INTC 
elseif  ind==2 
intplot2=odcorrha; 
610  save  INTM2  INTM 

save  INTK2  INTK 
save  DSTTC2  INTC 
elseif  ind==3 
iiitplot3=odcorrha; 
615  save  INTM3  INTM 

save  DSTTK3  INTK 
save  INTC3  INTC 
else 
iiUplot4=odcorrha; 
620  save  INTM4  INTM 

save  INTK4  INTK 
save  INTC4  INTC 
end 
end 
625 

if  meters 

\vaitbar_handle=waitbar(0,'Computing  Experimental  FRF'); 
else 
disp('Gettiiig  experimental  FRF') 
630        end 

HJEXP=[]; 

for  count=  1  :length(plotw) 
if  meters 
waitbar(count/length(plolw)); 
635  end 

z_exp=k_exp+j*plotw(count)*c_exp-plotw(count)A2*m_exp; 
h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
skyindex=l; 
640  for  indexl=l:aset_size 

for  index2=indexl  :aset_size 
red_holder(skyindex)=h_exp(  index  1  ,index2); 
skyindex=skyindex+ 1 ; 
end 
645  end 

H_EXP=[H_EXP  red_holder"]; 
end 

if  meters 
close(waitbar_handle) 
650        end 

figure(6); 

plot(plot\v/(2*pi),loglO(abs(H_EXP(ndx3d([l  aset_size*(aset_size+l )/2  length(plotw)], 

655  l,sk>Ted(cset_rel(mid_index),cset_rel(mid_index)),':')))),'r-',... 

plotvv/(2*pi),  Iogl0(abs(intplotl(ndx3d([aset_size  aset_size 
lengtli(plotw)],cset_rel(mid_mdex),cset_rel(mid_index),':')))),'g-',... 

plotw/(2*pi),  Iogl0(abs(intplot2(ndx3d([aset_size  aset_size 
lengtlnplotw)],cset_rel(mid_index),cset_rel(mid_mdex),':')))),'b:',... 
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plot\v/(2*pi),  Iogl0(abs(intplot3(ndx3d([aset_size  aset_size 
lengtli(plotw)],cset_rel(niid_index),cset_rel(niid_index),':')))),'k-.',... 

plotw/(2*pi),  Iogl0(abs(intplot4(ndx3d([aset_size  asetsize 
lengtli(plot\v)],cset_rel(mid_iudex),cset_rel(nud_iiidex),':')))),'m:') 
ylabeKf'HCintlstitcsetOnid^dexJXV^ntlstitcsetCmid^idex)),')  hi/lbf (loglO  of)']) 

xlabel('Omega  (Hz)') 

grid  on 

%tiUeCEXPERMENTAL  FRF  AND  CORRECTED  FRFS  USING  MODE  1  INT  SOLUTIONS') 

lih=sIegendCExp  FRF, ... 

[num2str(couiUctrl(l)),'  Pts  Corrected  FRF'], ... 

[nuni2str(countctrl(2)),'  Pts  Corrected  FRF'], ... 

[nuni2str(countctrl(3)),'  Pts  Corrected  FRF], ... 

jnumlstrCcountctrl^)),'  Pts  Corrected  FRF]); 
axes(hh); 

print  -dinfile  fig5_6 
figure(7); 
plot(plorw/(2*pi),loglO(abs(H_EXP(ndx3d([l  aset_size*(aset_size+l )/2  length(plolw)],... 

l,skyred(cset_rel(rriid_index),cset_rel(rnid_index)),':')))),'r-',... 

plotw/(2*pi),  Iogl0(abs(plotl(ndx3d([aset_sizeaset_size  lengtIi(plotw)],cset_rel(mid_index),cset_rel(mid_index),':')))),'g- 

plotw/(2*pi),  Iogl0(abs(iiitplotl(ndx3d([aset_size  aset_size 
lengtIi(plotw)],cset_rel(mid_mdex),cset_rel(nud_index),':')))),'b:') 
yIabel(['H(,,int2str(cset(mid_index)),',,,mt2str(cset(nud_index)),,)  in/lbf  (loglO  of)']) 

xlabel('Oraega  (Hz)') 
if  titles 

%  titleCEXPERMENTAL  FRF  AND  CORRECTED  FRFS  USING  MODE  1  OD&INTGRL  SOLTNS') 
end 
grid  on 
Wi=IegendCExp  FRF, ... 

[nuin2str(lenctrl(l)),'  Hz  Bandwidth  MAT  Corrected  FRF"], ... 

[num2str(leiictrl(  1 )),'  Hz  Bandwidtli  INT  Corrected  FRF]); 
axes(hh) 
if  pswitch==y 

print  -dcdj  color 
end 
print  -dinfile  fig5_7 

figured); 

plot(plotw/(2*pi),loglO(abs(H_EXP(ndx3d([l  aset_size*(aset_size+l )/2  length(plotw)],... 

1  ,skyred(cset_rel(mid_index),cset_rel(imd_index)),':')))),'r-',. . . 

plorw/(2*pi),  Iogl0(abs(plot2(ndx3d([aset_size  aset_size  lengtli(plotw)],cset_rel(mid_index),cset_rel(mid_index),':')))),'g- 

plohv/(2*pi),  Iogl0(abs(intplot2(ndx3d([aset_size  aset_size 
lengtli(plotw)],cset_rel(mid_index),cset_rel(mid_index),':')))),'b:') 
ylabel(  rH(',inastr(cset(nud_index)),',,,int2str(cset(nud_index)),')  in/lbf  (log  1 0  of)'] ) 

xlabel('Omega  (Hz)') 
if  titles 

%title(EXPERMENTAL  FRF  AND  CORRECTED  FRFS  USING  MODE  1  OD&INTGRL  SOLTNS') 
end 
grid  on 
hh=legendOExp  FRF, ... 

[nuni2str(lenctrl(2)),'  Hz  Bandwidth  MAT  Corrected  FRF'], ... 
[num2str(lenctrl(2)),'  Hz  Bandwidtli  INT  Corrected  FRF']); 
axes(hh) 
if  ps\\itch==y 
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720  print  -dcdj  color 

end 
print  -dinfile  fig5_8 

figure(9); 
725        plol(plotw/(2*pi),loglO(abs(H_EXP(ndx3d([l  aset_size*(aset_size+l)/2  length(plotw)],... 

1  ,skyred(cset_rel(mid_index),cset  jrel(mid_index)),':')))),'r-',. . . 

plotw/(2*pi),  Iogl0(abs(plot3(ndx3d([aset_sizeaset_size  length(plotvv)],cset_rel(mid_mdex),cset_rel(mid_index),':')))),'g- 

730  plotw/(2*pi),  Iogl0(abs(intplot3(ndx3d([aset_size  aset_size 

lengtli(plot\v)],cset_rel(mid_index),cset_rel(imd_index),':')))),'b:') 
ylabel(['H(',int2str(cset(nud_index)),,,',int2str(cset(mid_index)),')  in/lbf(loglO  of)']) 

xlabel('Omega  (Hz)') 
735        if  titles 

%title(*EXPERMENTAL  FRF  AND  CORRECTED  FRFS  USING  MODE  1  OD&INTGRL  SOLTNS") 
end 
grid  on 

Wi=legend(Exp  FRF', ... 
740  [num2str(lenctrl(3)),'  Hz  BandwidUi  MAT  Corrected  FRF'], ... 

[num2str(lenctrl(3)),'  Hz  Bandwidth  INT  Corrected  FRF]); 
&\es(hh) 
ifpswitch==y 
print  -dcdj  color 
745        end 

print  -dmfile  fig5_9 
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clear  , 

odivisions=3 

startloop^  1 

skiploop=  1 
5        endloop=4 

cd  incompl 

save  odcoiif  odi  visions  startloop  skiploop  endloop 

chap6_l 

chap6_2 
10        cd.. 

clear 

odivisions=l 

startloop=3 

skiploop=  1 
1 5        endloop=4 

cdcompl 

save  odconf  odi  visions  startloop  skiploop  endloop 

chapo_l 

chap6_3 
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20        %%%CHAP6_1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/0%%%%0/00/0%%% 
%%%%%%%%% 

%%%Decompose  DZ  into  DM,  DK,  &  IT/o%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%usmg  matrix  formulation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/0% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%%%0/o%%0/o 

25        %%%%%%%%%%%% 

clc 

clear 

closeall 

load  INT 
30        %ifpswitch=='y' 

whitebg('white') 

close 

%end 

h=0; 
35        fignum^l; 

use_antires=l; 

titles=0; 

mid_index=Tound(length(cset)/2); 

40        clear  L  Z_ANAL  cexp  conn  h_anal_red 

clear  kexto  mexp  mexto  temp_l_diags 

clear  truedamping  true_mass  true_stiffhess  z_anal  z_anal_red  z_exp 

clear  L  Z_ANAL  c_exp  conn  h_anal_red 

clear  kexta  kexto  mextO  temp_l_diags 
45 

%%%%%%%%%%%%%%%%%%%FORCETHECSET%%%%%0/o0/o0/o%0/o%%%%0/o0/o% 

%%%%%%%%%%%%SPATIALLYINCOMPLETE%0/o%%0/o%%%%%%%%%0/o0/o0/o0/o 

if  complete 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%0/o0/o%0/o0/o0/o0/o0/o% 
50        %%%%%%%%%%%%%%%SPATIALLYCOMPLETE%%0/o0/o%%%0/o0/o%0/o%%%0/o0/o0/o 
cset=[4  5  6  7  8  9  10  11]%%%%%%%%%%%%%%%%%%%%% 
cset_rel=cset%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
cset_size=lengdi(cset)%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
55        else 

%%%%%%%%%%%%SPATIALLYINCOMPLETE%%%%%%%%%%%%%%%%%% 
cset=[l  3  5  79  11  13  15]%%%%%%%%%%%%%%%%%%%%%% 
cset_rel=[l  2  3  4  5  6  7  8]%%%%%%%%%%%%%%%%%%%%% 
cset_size=lengtli(cset)%%%%%%%%%%%%%%%%%%%%%%%%% 
60        end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

load  exp 
65        odivisions=3 ; 

load  odconf 

nummodes=8    ; 

odlengtli=2.001*pi; 

firsttime=l; 
70        dw=odlengtli/(odivisions+l); 

dovv=(  freqtop-freqbottom)/(  fmeness- 1 ); 

owref=freqbottom:dow:freqtop; 

if  firsttime 
75  if  meters 

waitbar_handle=waitbar(0,'Computing  Experimental  FRF'); 
else 

disp('Getting  experimental  FRF') 
end 
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80  H_T_EXP=[]; 

for  count=l  :length(owref) 
if  meters 

waitbar(couiit/lengtli(owref)); 
end 
85  z_exp=k_exp+j  *owref(couiit)*c_exp-owref(couiit)A2  *m_exp; 

h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
%    skyiiidex-1; 
%    for  index  1 = 1  :aset_size 
90  %      forindex2=indexl:aset_size 

%         red_holder(  skyindex)=h_exp(  index  1 ,  index2 ); 
%         skyindex=skyindex+l; 
%      end 
%   end 
95  %   H_T_EXP=[H_T_EXP  redjiolder']; 

H_T_EXP=[H_T_EXP  h_exp(:)J; 
end 

if  meters 
close(waitbar_handle) 
100  end 

save_H_T_EXP  H_T_EXP 
clear  H_T_EXP 
%load  H_ANAL_RED; 
%uncorrha=H_ANAL_RED; 
1 05  %clear  H_ANAL_RED 

uncorrha=[]; 
if  meters 

waitbar_handle=waitbar(0,'Computing  Uncorrected  FRF'); 
else 
1 1 0  dispCGetting  Uncorrected  FRF') 

end 

load  stat 

for  i=l  :lengdi(owref) 
if  meters 
115  waitbar(i/lengtl\(o\vref)); 

end 

tempza=kstat+j*owref(i)*cstat-owref(i)A2*mstat; 
tempha=inv(tempza); 
uncorrha=[uncon"ha  tempha(:)]; 
120  end 

if  meters 

close(waitbarhandle) 
end 

save  uncorrha  uncorrha 
1 25  clear  uncorrha 

end 

load  H_T_EXP 
temp=[]; 

for  index=  1  :nummodes 
1 30  for  coord=  1  :aset_size 

if  index  ==  1 

antiresfreqs=fnid(owref>0  &  owref<omegax(index)); 
else 
antiresfreqs=find(owref>omegax(index-l)  &  o\vref<omegax(  index)); 
135  end 

antires=min(logl0(abs(H_T_EXP(ndx3d([aset_sizeaset_size  lengtli(ovvref)],. 

coord,coord,antiresfreqs))))); 
whre=fmd(logl0(abs(H_T_EXP(ndx3d([aset_sizeasel_size  lengU\(o\\Tef)],... 
coord,coord,antiresfreqs))))==antires), 
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140  lemp=owref(antiresfreqs) ; 

antiresfreq(coord,index)=temp(whre(  1 )); 
end 
end 

lemp=[]; 
145        for  index=l:nummodes 

resfreqs=find(owref>omegax(index)-dow  &  owref<omegax(indexyi-dow); 
res=max(log!0(abs(H_T_EXP(ndx3d([aset_sizeaset_size  length(owref)],... 

cset_rel(  1  ),cset_rel(  1  ),resfreqs))))); 
whre=fmd(logl0(abs(H_T_EXP(ndx3d([aset_sizeaset_size  length(owref)], 
1 50  cset_rel(  1  ),cset_rel(  1  ),resfreqs))))==res); 

temp=owref(resfreqs) ; 
resfreq(  index  )=temp(whre(  1 )); 
end 
clear  H  T  EXP 
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%clear  omegax 
%omegax=resfreq; 


for  loopindex=startloop:skiploop:endloop 
160  ODZ=[]; 

temp=[]; 
rdiag=[]; 
ow=[]; 

lo\vwer=o\vref(  1  )-.5*odlength; 
165  uppper=owref(l)+.5*odlength ; 

%  if  odivisions==  1 
%      o\v=[ow  owref(  1 )] 
%  else 

%      o\v=[ow  lowwer+dw:dw:uppper-dw]; 
170         %  end 

%  rdiag=[rdiag  ones(  1  ,odivisions*cset_size). *owref(  1 )]; 


for  index=  1  :loopindex 
175        %       lowwer=omega\(  index  )-.5*odlenglh ; 
%       uppper=omegax(index)+.5*odlength ; 

lo\v\ver=resfreq(index>.5*odlength ; 

uppper=resfreq(index)+.5*odlengtli ; 

1 80  o\v=[ow  lowwer+dw:dw:uppper-dw], 

rdiag=[rdiag  ones(l  ,odivisions*cset_size).*resfreq(index)]; 
if  use_antires 

for  coord=  1  :aset_size 
1 85  lowwer=antiresfreq(coord,index)-.5*odlengtli; 

uppper=antiresfreq(coord,index)+-.  5  *odlength; 
size(lowwer+d\v:dw:uppper-dw) 
ow=[ow  lowwer+dw:dw:uppper-dw] 
rdiag=[rdiag  ones(  1  ,odivisions*cset_size).*resfreq(index)]; 
190  end 

end 
end 

R=diag(rdiag); 
load  stat 
195  if  meters 

\vaitbar_handle=waitbar(0,'ComputmgDK/DM/DC); 
else 

dispCGetting  DK/DK4/DC) 
end 
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200 

for  count=  1  :length(ow) 
if  meters 

waitbar(count/length(ow)); 
end 
205  z_anal_red=kstat-ow(count)A2  *mstat; 

h_aiial_red=iiiv(z_anal_red); 

z_exp=k_exp+j  *ow(count)*c_exp-ow(  count  )A2  *m_exp; 
h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
210  hacc=h_anal_red(cset_rel,cset_rel); 

hxcc=h_exp(cset_rel,cset_rel); 
dz=inv(inv(inv(hacc)*(hacc-lixcc)*inv(hacc))-hacc); 
ODZ=[ODZ;  dz]; 

temp=[temp;  eye(cset_size)  -ow(count)A2*eye(cset_size)  +j*ow(count)*eye(cset_size)]; 
215  end 

if  meters 

close(waitbar_handle) 
end 

DKDMDC 1 =inv(temp'*inv(R)*temp)*teinp'*inv(R)*ODZ; 
220  %DKDMDC  1 =temp\ODZ ; 

ODKl=DKDMDCl(l:cset_size,:) ; 

ODM 1  =DKDMDC  1  (cset_size+ 1 :2  *cset_size, : ) ; 

ODCl=DKDMDCl(2*cset_size+l  :3*cset_size,:) ; 

clear  hacc  lixcc  h_exp  zexp  z_anal_red  h_anal_red  dz 
clear  temp  ODZ 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o0/o%%%% 

%%%%%%%%%%%% 

load  star, 
odcorrha=[]; 

kstat(cset_rel,cset_rel)=kstat(cset_rel,cset_rel)+ODKl; 
mstat(cset_rel,cset_rel)=mstat(cset_rel,cset_rel)+ODM  1 ;; 
cstat(cset_rel,cset_rel)=cstat(cset_rel,cset_rel)+ODC  1 ; 

if  meters 

waitbar_handle=waitbar(0,'Computing  Corrected  FRF); 
else 

disp('Getting  DK/DM/DC') 
end 
for  i=l:length(owref) 

if  meters 
waitbar(i/length(owref)); 

end 

tempza=kstat+j  *owref(i)*cstat-owref(i)A2  'mstat; 

tempha=inv(tempza); 

odcorrha=[odcorrha  templia(:)]; 
end 
if  meters 

close(vvaitbar_liandle) 
end 
if  loopindex  ==  1 

odcorrha  1  =odcorrha; 

save  odl  odcorrha  1 

clear  odcorrhal 
elseif  loopindex  ==  2 

odcorrha2=odcorTha; 

save  od2  odcorrha2 
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260  clear  odcorrha2 

elseif  loopindex  ==  3 
odcorrha3=odcorrha; 
save  od3  odcorrha3 
clear  odcorrha3 
265  elseif  loopindex  ==  4 

odcorrha4=odcorrha; 
save  od4  odcorrha4 
clear  odcorrha4 
elseif  loopindex  ==  5 
270  odcorrha5=odcorrha; 

save  od5  odcorrha5 
clear  odcorrha5 
end 

275  clear  tempha  tempza  kcorrected  mcorrected  ccorrected 

clear  ODM1  ODC1  DKDMDC1  ODZ  temp  kstat  mstat 

load  uncorrha 
load  H_T_EXP 
280  if  complete 

plotwluchcoord=[cset(mid_index)]; 
else 

plotwhichcoord=[cset(mid_index)]; 
end 
285  for  plotindex=  1  :length(plorwhichcoord) 

which_coord=plotwhichcoord(plotindex); 
h=figure(h+l); 
fignum=fignum+l  ; 
subplot(2,l,l); 
290  plot(owref/(2*pi),... 

loglO(abs(H_T_EXP(ndx3d([aset_size  aset_size  lengtli(ovvref)],... 
which_coord,which_coord,':,)))),'r-',. . . 
owref/(2*pi),... 
Iogl0(abs(odcorrha(ndx3d([aset_size  asetsize  length(owref)],... 
295  wluch_coord,which_coord,':')))),'g— ',. . . 

o\vref/(2*pi),... 
Iogl0(abs(uncorrha(ndx3d([aset_size  aset_size  length(o\vref)],-- 
wlhch_coord,\vluch_coord,':')))),'b-.') 
ylabel(['H(',int2str(aset(which_coord)),... 
300  7,int2str(aset(which_coord)),')  in/lbf  (loglO  of)']) 

hold  on 
v=axis; 

%xlabel('Omega  (Hz)') 
if  titles 
305  if  use_antires 

title([int2str(odivisions),'  PT  MATRIX  SOLTNS  WTTH  ANTIRESONANCES  WEIGHTED  LEFf  TO  RIGHT]) 
else 

titleamastrCodivisions),'  PT  MATRIX  SOLTNS  WEIGHTED  LEFT  TO  RIGHT]) 
end 
310  end 

plot(resfreq(  1  :loopindex)/(2*pi),... 
ones(  1  ,loopindex)*(  v(3)+abs(v(4)-v(3))*.98),'k+') 
%       ifuse_antires 

%      %plot(diag(mitiresfreq(l:lc)opmdex,l:loopmdex))/(2*pi),ones(lJoopmde%x)*(v(3)+abs(v(4)-v(3))*.98),'k*') 
315        %       end 
grid  on 

hh=legendCExperimental  FRF','Corrected  MAT  Anal  FRF/Uncorrected  Anal  FRF,1ncluded  Modes'); 
axes(lili) 
hold  off 
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320  %ifpswitch=='y' 

%    prtfigl(fignum) 
%    delete(h) 
%      %  else 
%  delete(h) 

325  %  end 

subplot(2,l,2) 

tiietop=find(abs(owref-ow(length(ow))*  1 . 5  *ones(  1  ,length(owref)))==niin(abs(o\vref- 
ow(length(ow))*  1 .5*ones(  1  ,length(owref))))); 

330  plol(owref(l  :thetop)/(2*pi),... 

loglO(abs(H_T_EXP(ndx3d([aset_size  aset_size  lenglh(ovvref)],... 
which_coord,wluch_coord,  1  :thetop)))),'r-',. .  • 
owref(l  :thetop)/(2*pi),... 
Iogl0(abs(odcorrha(ndx3d([aset_size  aset_size  length(owref)],... 
335  wluch_coord,which_coord,  1  :thetop)))),'g— ',... 

owref(l:thetop)/(2*pi),... 
Iogl0(abs(uncorrha(iidx3d([aset_size  aset_size  length(owref)],... 
wluch_coord,wluch_coord,  1  :thetop)))),t>-.') 
ylabel(['H(,,int2str(aset(wluch_coord)),',*,int2str(aset(wliich_coord)),')  iii/lbf  (log  1 0  of)']) 
340  hold  on 

v=axis; 

xlabeK'Omega  (Hz)*) 
if  titles 
if  use_antires 
345  title([int2str(odi visions),'  PT  MATRIX  SOLTNS  WITH  ANTIRESONANCES  WEIGHTED  LEFT  TO  RIGHT]) 

else 

tideaint2str(odivisions),'  PT  MATRIX  SOLTNS  WEIGHTED  LEFT  TO  RIGHT]) 
end 
end 
350  plot(resfreq(  1  :loopindex)/(2*pi),ones(  1  ,loopindex)*(v(3)+abs(v(4)-v(3))*.98),'k+') 

%       ifuse_antires 

%    %plot(antiresfreq(which_coord,  1  :loopindex)/(2*pi),ones(  1  ,loopindex)*(v(%3>i-abs(v(4)-v(3))*.98),,k*') 
%       end 
grid  on 
355  %hh=legendCExperimental  FRF','Corrected  MAT  Anal  FRF','Uncorrected  Anal  FRF^Tncluded  Modes'); 

%axes(hh) 
hold  off 

if  loopmdex  —  2 
print  -dcdjcolor 
360  print  -dinfile  fig6_l 

elseif  loopmdex  ==  3 
ifodivisions==  1 
print  -dcdjcolor 
if  complete 
365  print  -dmfile  fig6_9 

else 

print  -dmfile  fig6_6 
end 
end 
370  elseif  loopindex  ==  4 

print  -dcdjcolor 
print  -dmfile  fig6_2 
end 
%     if  pswitch==y 
375  %        prtfigl(fignum) 

%        delete(h) 
%    else 
%        delete(h) 
%    end 
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380  end      ' 

delete(h) 
clear  H_T_EXP 
clear  odcorrha 
clear  uncorrha 

385        end 
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%%%%CHAP6J.M%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/0%%%%%%%% 
%%%%%%%%% 

%%%  Decompose  DZ  into  DM,  DK,  &  DC  %%%Vo%%%%%%%%%%%%%y0%%Vo%%%%%%%%VoVo%%%%%% 
390        %%%  using  matrix  formulation       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%o/0%o/0o/0o/oo/o 

mm%%mm%%%m%%%%m%%%m%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%% 

clc 

clear 
395        closeall 

load  INT 

whitebg(\vhite') 

close 

h=0; 
400        fignum=l; 

use_antires=l 

titles=0 

mid_index=round(length(cset)/2); 

405       clear  L  Z_ANAL  c_exp  conn  h_anal_red 

clear  kexto  mexp  mexto  temp_l_diags 

clear  true_damping  true_mass  true_stiffness  z_anal  z_anal_red  z_exp 

clear  L  Z_ANAL  c_exp  comi  h_anal_red 

clear  kexta  kexto  mextO  temp_l_diags 
410 

%%%%%%%%%%%%%%%%%%%FORCETHECSET%%%%%0/o0/o0/o0/o0/o0/o0/o%0/o%%% 

%%%%%%%%%%%%SPATIALLYINCOMPLETE%%0/o%0/o%%%0/o0/o0/o0/o0/o%0/o0/o%0/o 

if  complete 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o%%%%%0/o%0/o%%0/o%%% 

4 1 5        %%%%%%%%%%%%%%%SPATIALL Y  COMPLETE%%%%%%%%%%%%%%%%% 
cset=[4  5  6  7  8  9  10  1 1]%%%%%%%%%%%%%%%%%%%%% 
csetj-el=K:set%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
cset_size=length(cset)%%%%%0/o%0/o%0/o%%0/o%0/o%%0/o0/o0/o0/o0/o%%0/o0/o% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%0/o%0/o%%0/o%%0/o0/o%%%% 

420        else 

%%%%%%%%%%%%SPATIALLYINCOMPLETE%0/o0/o0/o%0/o%%%0/o0/o%0/o0/o0/o%%0/o 

cset=[l  3  5  79  11  13  15]%%%%%%%%%%%%%%%%%%%%%% 

cset_rel=[l  2  3  4  5  6  7  8]%%%%%%%%%%%%%%%%%%%%% 

cset_size=lengm(cset)%%%%%%%%%%%%%%%%%%%%%%%%% 
425        end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%SPATIALLYCOMPLETE%%%%%%%%%%%%%%%%% 
430        %cset=[7  8  9  10  11  12  13  14]%%%%%%%%%%%%%%%%%%%%% 

%cset_rel=cset%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%cset_size=lengUi(cset)%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

435        load  exp 

odivisions=3  ; 

load  odconf 

nununodes=9     ; 

odlength=2*pi ; 
440        firsttime=l; 

dw=odlengdi/(odivisions-t-l ); 

dovv=(freqtop-freqbottom)/(  fineness- 1 ); 

owref=freqbottom:dow:freqlop; 

445        iffirsttime 
if  meters 
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waitbar_handle=waitbar(0,'Computing  Experimental  FRF'); 
else 
disp('Gettiiig  experimental  FRF') 
450  end 

H_T_EXP=[]; 
for  count=  1  :length(owref) 
if  meters 
455  waitbar(count/length(ovvref)); 

end 

z_exp=k_exp+j*owref(count)*c_exp-o\vref(  count  )A2*m_exp; 
h_exp=inv(z_exp); 
h_exp=h_exp(aset,aset); 
460        %      skyindex=l; 

%      for  indexl = 1  :aset_size 
%         for  mdex2=indexl:aset_size 
%  red_holder(skyindex)=h_exp(  index  1  ,index2); 

%  skyindex=skyindex+l; 

465        %         end 
%     end 

%     H_T_EXP=[H_T_EXP  redjiolder']; 
H_T_EXP=[H_T_EXP  h_exp(:)]; 
end 


470 


475 


if  meters 

close(waitbar_handle) 
end 


save  H_T_EXP  H_T_EXP 
clear  H_T_EXP 
end 
480        %load  H_ANAL_RED; 

%uncorrha=H_ANAL_RED; 
%clear  H_ANAL_RED 
if  firsttime 
uncorrha=[]; 
485  if  meters 

waitbar_handle=waitbar(0,'Computing  Uncorrected  FRF'); 
else 

disp('Getting  Uncorrected  FRF') 
end 
490 

load  stat 

for  i=  1  :length(owref) 
if  meters 
waitbar(i/lengtli(ovvTeO); 
495  end 

tempza=kstat+j  *owref(i)*cstat-owref(i)A2  *mstat; 
tempha=iiiv(tempza); 
uncorrha=[uncorrha  tempha(:)]; 
end 
500 

if  meters 

close(wailbar_handle) 
end 

505  save  uncorrha  uucorrha 

clear  uncorrha 
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end 

load  H_T_EXP 
510        for  iiidex=  1  :nummodes 
for  coord=  1  :aset_size 
if  index  —  1 

antiresfreqs=find(owref>0  &  owref<omegax(index)); 
else 
515  anthesfreqs=find(owref>omegax(index-l )  &  owref<omegax(index)); 

end 
aiitires=inin(logl0(abs(H_T_EXP(ndx3d([l  aset_size*(aset_size+l )/2  length(owref)],. 

l,skyred(coord,coord),aiitiresfreqs))))); 
wlire=fmd(loglO(abs(H_T_EXP(ndx3d([l  aset_size*(aset_size+l  )/2  length(owref)],... 
520  1  ,skyred(coord,coord),antiresfreqs))))==aiitires); 

lemp=owref(antiresfreqs) ; 
antiresfreq(coord,index)=temp(whre); 
end 
end 
525 

for  index-  lmummodes 

resfreqs=find(owref>omegax(index)-dow  &  owref<omegax(index)+-dow); 
res=max(logl0(abs(H_T_EXP(ndx3d([aset_sizeaset_size  length(owref)],... 
cset_rel(  1  ),cset_rel(  1  ),resfreqs))))); 
530  wlire=find(loglO(abs(H_T_EXP(ndx3d([aset_size  aset_size  length(owref)],... 

cset_rel(  1  ),cset_rel(  1  ),resfreqs))))==res); 
temp=o\vref(resfreqs) ; 
resfreq(  index  )=temp(whre); 
end 
535 

clear  H_T_EXP 
%clear  omegax 
%omegax~resfreq; 

540        for  loopindex=startloop:skiploop:endloop 

subdivisions=9;  %  should  be  odd  for  sinipson's  rule 

intlengtli=2*pi;  %  1  Hz  bandwidth 

dW=intlength/subdivisions;  %  sampling  freq  =.25  HZ 

W=0; 

545  Wl=[]; 

vveiglit=[]; 

low\ver=o\vref(  1  )-.5*odlength; 

uppper=owref(  1  )+.5*odlength; 

W=[W  lo\vwer+dW:dW:uppper-dW]; 
550  W1=[W1  low\ver+dW:dW:uppper-dW]; 

weigh  t=  [weight  ones(size(lowwer+dW:dW:uppper-dW))]; 

for  mdex=  1  :loopindex 
lowwer=omegax(index)-.5*intlengtli; 
555  uppper=omegax(  index  )+.  5*intlengtli; 

W=[W  lowwer+dW:dW:uppper-dW]; 
W1=[W1  lowwer+dW:dW:uppper-dW]; 
weight=  [weiglit  ones( size( lowwer+dW:dW : uppper-d W))] ; 
lo\ver=antiresfreq(index)-.5*intlengtli; 
560  upper=antiresfreq(  index  }+.5  *mtlength; 

W=[W  lowwer+dW:dW:uppper-dW]; 
W1=[W1  lowwer+dW:dW:uppper-dW]; 
weight=[weight  oues(size(lowwer+dW:dW:uppper-dW))]; 
end 


565 


%%computemtegrals%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o0/o0/o%%0/o%0/o%%%%0/o 
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DZ_R=[]; 

pz_i=D; 

if  meters 
570  waitbar_haiidle=waitbar(0,'Computing  DZ'); 

else 

disp('Computing  DZ'); 
end 

load  stat 
575  for  count=  1  :lengdi(  W) 

if  meters 

waitbar(count/length(W)); 
end 

z_anal_red=kstat/(j  *  W(count)>i-cstat/(j  *  W(  count  ))+-j  *  W(coimt)*mslat; 
580  h_anal_red=inv(z_anal_red); 

hacc=h_anal_red(cset_rel,cset_rel); 

z_exp=k_exp/(j  *  W(count))+c_exp/(j  *  W(count))+j  *  W(count)*m_exp; 
h_exp=inv(z_exp); 
585  h_exp=h_exp(aset,aset); 

lrxcc=h_exp(cset_rel,cset_rel); 

dz=inv(mv(inv(hacc)*(hacc-lixcc)*inv(hacc))-hacc); 
dz_r=real(dz); 
590  dz_i=imag(dz); 

DZ_R=[DZ_R  dz_r(:)]; 
DZ_I=pZ_I  dz_i(:)]; 
end 

595  if  meters 

close(waitbar_handle) 

end 

W1=1./W1; 

[INTK,  INTM,  INTC]=intsub(DZ_LDZ_R,W,Wl,weiglit,cset_size); 
600  intcorrha=[]; 

clear  hacc  hxcc  h_exp  zexp  z_anal_red  hanalred  dz 

clear  temp  ODZ 

605        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o%%%% 
%%%%%%%%%%%% 
load  stat; 
odcorrha=[]; 

610  kstat(cset_rel,cset_rel)=kstat(cset_rel,cset_rel>i-INTK; 

mstat(cset_rel,cset_rel)=mstat(cset_rel,cset_rel>+-INTM; 
cstat(cset_rel,cset_rel)=cstat(cset_rel,cset_rel)+INTC; 

if  meters 
615  \vaitbar_handle=waitbar(0,'Computing  Corrected  FRF'); 

else 

disp('Getting  DK/DM/DC') 
end 

620  for  i=  1  :lengtli(owref) 

if  meters 

waitbar(i/length(ovvref)); 
end 

tempza=kstat+j*owTef(i)*cstat-owref(i)A2*mstat; 
625  tempha=inv(tenipza); 

odcorrha^odcoiTha  tcmpha(:)]; 
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end 

if  meters 

close(waitbar_handle) 
end 

if  loopindex  ==  1 

lodcorrha  1  =odcorrha; 

save  lodl  lodcorrha  1 

clear  lodcorrha  1 
elseif  loopindex  =  2 

lodcorrha2=odcorrha; 

save  lod2  lodcorrha2 

clear  Iodcorrha2 
elseif  loopindex  ==  3 

lodcorTha3=odcorrha; 

save  lod3  lodcorrha3 

clear  lodcorrha3 
elseif  loopindex  ==  4 

lodcorrha4=odcorrha; 

save  lod4  lodcorrha4 

clear  lodcorrha4 
elseif  loopindex  ==  5 

lodcorrha5=odcorrha; 

save  lod5  lodcorrha5 

clear  lodcorrha5 
end 


clear  tempha  tempza  kcorrected  mcorrected  ccorrected 
clear  INTM INTC  INTK  temp  kstat  mstat  cstat 
load  uncorrha 
load  H_T_EXP 

if  complete 

plot\vliichcoord=[cset(nud_index)]; 
else 

plotwhichcoord=[cset(mid_index)]; 
end 

for  plotindex=  1  :length(plotwhichcoord) 
wluch_coord=plotwhichcoord(plotindex); 
h=figure(h+l); 
subplot(2,l,l) 
plol(owref/(2*pi),... 

loglO(abs(H_T_EXP(ndx3d([aset_size  aset_size  lengtli(owref)], 

wluch_coord,which_coord,':')))),'r-',. . . 
o\vref/(2*pi),... 

Iogl0(abs(odcorrha(ndx3d([aset_size  aset_size  lengtli(owref)],... 

wluch_coord,which_coord,':,)))),'g-',.. 
owref7(2*pi),... 

Iogl0(abs(uncorrha(ndx3d([aset_size  aset_size  lengtli(o\AreO], •■• 

wluch_coord,wluch_coord,':')))),'b-.') 
ylabel(['H(,,int2str(aset(wluch_coord)),. . . 

7,int2str(aset(wluch_coord))/)  in/lbf  (log  10  o0']) 
hold  on 
v=axis; 

%xlabel('Omega  (Hz)') 
if  titles 

if  use  antires  ==  'on ' 
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title([int2str(odi visions),'  PT  INTEGRAL  SOLTNS  WITH  ANTIRESONANCES  WEIGHTED  LEFT  TO  RIGHT]) 
else 
title([int2str(odivisions),'  PT  INTEGRAL  SOLTNS  WEIGHTED  LEFT  TO  RIGHT]) 
690  end 

end 

plot(omegax(  1  :loopindex)/(2*pi),. . . 
ones(  1  ,loopindex)*(v(3)+abs(v(4>v(3))*.98),'k+') 
%      if  use_antires 

695        %         %plot(aiitiresfreq(wliich_coord,l:loopindex)/(2*pi),ones(l,loopindex)*(v(%3)+abs(v(4>v(3))*.98),'k*') 
%       end 
grid  on 

hh=legendCExperimental  FRF/Corrected  INT  Anal  FRF^TJncorrected  Anal  FRF /Included  Modes'); 
axes(hh) 
700  hold  off 

%ifpswitch==y 
%   prtfigl(fignuin) 
%    delete(h) 
%else 
705  %   delete(h) 

%end 

subplot(2,l,2) 

thetop=find(abs(owref-W(length(W))*  1 .3*ones(  1  ,length(owref)))==rnin(abs(o\vref- 
W(lengtIi(W))*  1 .3*ones(  1  ,length(owref))))); 
710  %h=figure(h+l); 

%  fignum=fignum+l  ; 
plot(owref(  1  :thetop)/(2*pi),... 
logl0(abs(H_T_EXP(ndx3d([aset_size  asetsize  length(owref)],... 
which_coord,which_coord,  1  :thetop)))),'r-',. . . 
7 1 5  owref(  1  :thetop)/(2*pi),. . . 

Iogl0(abs(odcorrha(ndx3d([aset_size  aset_size  length(o\vref)],.- 
which_coord,  which_coord,  1 :  thetop )  ))),'g-', . . . 
owref(  1  :thetop)/(2*pi),... 
logl0(abs(uncorTha(ndx3d([aset_size  asetsize  Iengtli(owre0], ... 
720  wliich_coord,wliich_coord,  1  itlietop))))^.') 

ylabel(['H(',int2str(aset(wluch_coord)),',,,int2str(aset(\vhich_coord));)  ui/lbf  (loglO  ol")']) 
hold  on 
v=axis; 

xlabel('Omega  (Hz)') 
725  if  titles 

if  use_antires  ==  'on ' 

title([mt2str(odivisions)/  PT  INTEGRAL  SOLTNS  WITH  ANTIRESONANCES  WEIGHTED  LEFT  TO  RIGHT]) 
else 
title([int2str(odi visions),'  PT  INTEGRAL  SOLTNS  WEIGHTED  LEFT  TO  RIGHT]) 
730  end 

end 

plot(omegax(  1  :loopindex)/(2  *pi),ones(  1  ,loopmdex)*(  v(3  )+abs(  v(4  >v(  3  ))* .  98),"k+') 
%      if  use_antires 

%  %plot(antiresfreq(which_coord,l:loopindex)/(2*pi),ones(l,loopindex)*(v(%3)+abs(v(4)-v(3))*.98),'k*') 

735        %       end 
grid  on 

%hh=legendCExperimental  FRF','Corrected  INT  Anal  %RF','Uncorrected  Anal  FRF',Tncluded  Modes'); 
%axes(lili) 

740  if  loopindex  =  2 

print  -dcdj  color 

print  -dinfile  fig6_3 

h=figure(h+l); 

hold  off 
745  load  od2 

plot(owTef(  1  :thetop)/(2*pi),... 
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loglO(abs(H_T_EXP(ndx3d([aset_size  aset_size  length(owref)],- 
\vliich_coord,wliich_coord,l:tlietop)))),'r-',..- 
owref(  1  :thetop)/(2*pi),... 
750  Iogl0(abs(odcorrha(udx3d([aset_size  aset_size  lengUi(owref)],-.. 

wliich_coord,vvliich_coord,l:tlietop)))),'g-',... 
o  wref(  1  :tlietop)/(2  *pi),. . . 
loglO(abs(odcorTha2(ndx3d([aset_size  aset_size  length(owref)],-- 
wluch_coord,wliich_coord,  1  itlietop)))),!)-.',. . . 
755  owref(l:thetop)/(2*pi),... 

logl0(abs(uncoiTha(ndx3d([aset_size  aset_size  length(owref)],... 
which_coord,which_coord,  1  :thetop)))),'m:') 
ylabel([H(',iiit2str(aset(wluch_coord)),V,iiit2str(aset(wliich_coord)),')  in/lbf(loglO  of)']) 
hold  on 
760  v=axis; 

xlabel('Omega  (Hz)') 
if  titles 
if  use_antires  ==  'on ' 
title([int2str(odivisions),'  PT  INTEGRAL  SOLTNS  WITH  ANTIRESONANCES  WEIGHTED  LEFT  TO  RIGHT]) 
765  else 

titleaii^strCodivisions),*  PT  INTEGRAL  SOLTNS  WEIGHTED  LEFT  TO  RIGHT]) 
end 
end 

plot(omegax(l:loopindex)/(2*pi),ones(l,loopindex)*(v(3)+abs{v(4)-v(3))*.98),'k+') 
770        %         ifuse_antires 

%  %plot(antiresfreq(which_coord,  1  :loopindex)/(2*pi),ones(  1  ,loopindex)*(v(%3>fabs(v(4>v(3))*.98),'k*') 

%         end 
grid  on 

hh=legend('Experimental  FRF'/Corrected  INT  Anal  FRF,'Corrected  MAT  Anal  FRF',TJncorrected  Anal  FRF','Includcd 
775        Modes'); 

axes(hh) 
hold  off 
print  -dcdj  color 
print  -dinfile  fig6_5 
780  elseif  loopindex  =  3 

if  odivisions==  1 
print  -dcdjcolor 
print  -dinfile  fig6_7 
h=figure(h+l); 
785  hold  off 

loadod3 

plot(owref(  1  :tiietop)/(2  *pi),. .  . 
loglO(abs(H_T_EXP(ndx3d([aset_size  aset_size  length(owret)],... 
wluch_coord,wlnch_coord,  1  :thetop)))),'r-',. . . 
790  owrefU:Uietop)/(2*pi),... 

Iogl0(abs(odcorrha(ndx3d([aset_size  aset_size  lengtli(o\\Te0],... 
\vliich_coord,wluch_coord,  1  :thetop)))),'g— ',... 
owref(  1  :thetop)/(2*pi),... 
Iogl0(abs(odcorrha3(ndx3d([aset_size  aset_size  lengtli(o\vre0],-- 
795  which_coord,wluch_coord,  1  itlaetop)))),!)-.',. . . 

owref(l:thetop)/(2*pi),... 
Iogl0(abs(uncorrha(ndx3d([aset_size  aset_size  Iengtli(owre0],... 
\vliich_coord,wluch_coord,  1  :thetop)))),'m:') 
ylabel(['H(',int2str(aset(wluch_coord)),V,int2str(aset(which_coord)),')  in/lbf  (loglOot)']) 
800  hold  on 

v=axis; 

xlabelCOmega  (Hz)') 
if  tides 
if  use_antires  ==  'on ' 
805  title([mt2str(odi visions),'  PT  INTEGRAL  SOLTNS  WITH  ANTTRESONANCES  WEIGHTED  LEFT  TO  RIGHT]) 

else 
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title([int2str(odi visions),'  PT  INTEGRAL  SOLTNS  WEIGHTED  LEFT  TO  RIGHT]) 
end 
end 
810  plot(omegax(  1  :Ioopindex)/(2*pi),ones(  1  ,loopindex)*(v(3)+abs(v(4)-v(3))*.98),,k+') 

%  if  use_antires 

%  %pIot(antiresfreq(\vliich_coord,l  :loopindex)/(2*pi),ones(  1  ,loopindex)*(v(%3)+abs(v(4>v(3))*.98),'k*') 

%  end 

grid  on 
8 1 5  hh=legendCExperimental  FRF/Corrected  INT  Anal  FRFVUncorrected  Anal  FRI^Tncluded  Modes'); 

axes(lih) 
hold  off 
print  -dcdj  color 
print  -dmfile  fig6_8 
820 

end 
elseif  loopindex  ==  4 
print  -dcdj  color 
print  -dmfile  fig6_4 
825  end 

%       ifpswitch=='y' 
%         prtfigl(fignuni) 
%         delete(h) 
%     else 
830         %        delete(h) 
%     end 
end 

delete(h) 
clear  H_T_EXP 
835  clear  odcorrha 

clear  uncorrha 
end 

loadH_T_EXP 
840  load  odl 
loadod2 
loadod3 
load  od4 
load  uncorrha 
845        %load  lod5 

h=figure(h+l); 
plol(owref/(2*pi),... 
loglO(abs(H_T_EXP(ndx3d([aset_size  aset_size  length(owref)]v 
vvhich_coord,wliich_coord,':')))),'r-',... 
850  o\\Tef/(2*pi),... 

Iogl0(abs(odcorrhal(ndx3d([aset_size  aset_size  lengtli(owre0],-- 
\vluch_coord,\vhich_coord,':')))),'g-',. . . 
owref/(2*pi),... 
Iogl0(abs(odcorrha2(ndx3d([aset_size  aset_size  length(ov\Tef)],-- 
855  which_coord,which_coord,':')))),'b-.',. . . 

o\vref/(2*pi),... 
loglO(abs(odcorTha3(ndx3d([aset_size  aset_size  length(owref)],--- 
wluch_coord,\vliich_coord,':')))),'m— ',. . . 
o\vref/(2*pi),... 
860  loglO(abs(odcorTha4(iidx3d([aset_size  aset_size  length(o\vref)],-- 

wluch_coord,\vluch_coord,':')))),'c-.',. . . 
o\vref/(2*pi),... 
Iogl0(abs(uncorrha(ndx3d([aset_size  aset_size  length(owref)],... 
whichcoord^vliichcoord,':')))),^-.') 
865  yiabeKPHCii^st^ase^whicl^coord)),... 

V,int2str(aset(\vluch_coord),),')  in/lbf  (loglO  of)']) 
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hold  on 
v=axis; 

%xlabel('Omega  (Hz)') 
870  if  titles 

if  useantires  ==  'on ' 

title([int2str(odivisions),'  PT  INTEGRAL  SOLTNS  WITH  ANTTRESONANCES  WEIGHTED  LEFT  TO  RIGHT]) 
else 
litle([int2str(odivisions),'  PT  INTEGRAL  SOLTNS  WEIGHTED  LEFT  TO  RIGHT]) 
875  end 

end 

%plot(oraegax(l:loopindex)/(2*pi),... 
%     ones(l,loopindex)*(v(3)+abs(v(4)-v(3))*.98),,k+') 
%  if  use_antires 
880  %     %plot(antiresfreq(which_coord,l:loopindex)/(2*pi),ones(l,loopindex)*(v(0/o3)+abs(v(4)-v(3))*.98),'k*') 

%       end 
v=axis; 

axis([10  300  v(3)  v(4)]); 
grid  on 
885  Ui=legend(Experimental  FRF','1  mode  MAT  solution','2  mode  MAT  soluUon','3  mode  MAT  solution','4  mode  MAT 

solutioi^TJncorrected  Anal  FRF); 
axes(hh) 
hold  off 

%if  pswitch==y 
890  %   prtfigl(fignum) 

%    delete(h) 
%else 

%    delete(h) 
%end 


895 


900 


print  -dcdj  color 
print  -dmfile  fig6_12 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%o/0o/0o/o%%%%%o/0%o/0%o/o 

%%%%%%%%%% 

%  PROGRAM:  SETUP.M  % 

%  INTTILIZES  DATA  FILE  FOR  A  FINTTE  ELEMENT  ANALYSIS  % 

5        %  VARIABLES  SAVED  TO  FILE  SETUP.MAT  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%o/0o/o%%%%%%%%o/0o/o%%%%%o/o%%%%o/o%o/o%%%% 

%%%%%%%%%% 

clc 
10         clear 

L=input(Enter  the  total  length  of  the  beam: '); 
numeI=inputCNum  of  elements: '); 
dof_node=iiiput('uum  dof  per  node: '); 
area=input('area  of  beam: '); 
1 5        eeii=inputCEI  of  beam: '); 

%ee=input('modulus  for  beam: '); 
pho=input('mass  density  for  beam: '); 
numdof=dof_node*(numel+ 1 ) 
clc 
20        while  1 

conforce=input(Enter  #  of  concentrated  loads: '); 
if  -isnan(conforce) 

break 
end 
25        end 

for  i=l:conforce 
forcepos(i)=input(['Beam  position  of  force  ',num2str(i),' : ']); 
forcesiz(i)=input([Magnitude  of  force  ',num2str(i),' : ']); 
end 
30        force=zeros(numdof,l); 
for  i=l:conforce 
temppos=ceil(forcepos(i)/(L/numel)); 

force(dof_node*temppos+l  )=force(dof_node*temppos+l  )+.5*forcesiz(i); 
force(dof_node*temppos+3)=force(dof_node*temppos-f-3>t-.5*forcesiz(i); 
35 

end 
while  1 

lmass=inputCEnter  #  of  Lumped  masses: '); 
if-isnan(lmass) 
40  break 

end 
end 

fori=l:lmass 
masspos(i)=mput(['Beam  position  of  Mass  ',num2str(i),' : ']); 
45  masssiz(i)=input([Magnitude  of  Mass  '4ium2str(i),' : ']); 

end 

lumpmass=zeros(numel,  1 ); 
for  i=l:lmass 
temppos=round(masspos(i)/(L/numel)); 
50  lumpmass(temppos)=lumpmass(temppos)+masssiz(i); 

end 

lumpmass' 
wlule  1 

lspring=mputCEnter  #  of  Lumped  Springs: '); 
55  if~isnan(lspring) 

break 
end 
end 

for  i=  1  :lspring 
60  springpos(i)=input([T3eam  position  of  Spring  ',num2str(i),' : ']), 
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springsiz(i)=iiiput(['Spring  constant  of  Spring  ',nuni2str(i),' : ']); 
end 

liunpspring=zeros(numel,l ); 
fori=l:lspring 
65  temppos=round(springpos(i)/(L/nuineI)); 

luinpspring(temppos)=luinpspring(tenippos)+springsiz(i); 
end 

lumpspring' 
conn=[l,2]; 
70        for  i=2:numel 

conn=[conn;i,i+l]; 
end 
conn 
while  1 
75  bc=['pinned-pinned     ' 

'clamp-clamp       ' 
'left  guided  clamp ' 
'right  guided  clamp' 
'cantilevered      ' 
80  'free-free  ']; 

clc 

help  bctext 

n=input('Select  a  boundary  condition: '); 
if((n>0)&(n<7)) 
85  break 

end 
end 

ifn==l 
bc=«pp'; 
90        elseifn==2 
bc='cc'; 
elseif  n  =  3 

bc=lc'; 
elseif  n  ==  4 
95  bc='rc'; 

elseif  n==5 

bc='cl' 
else 
bc='ff; 
100        end 

clear  i 

clear  temppos 

clear  conferee 
1 05        clear  forcepos 

clear  forcesiz 

clear  lmass 

clear  inasspos 

clear  masssiz 
110        clear  lspring 

clear  springpos 

clear  springsiz 

save  setup.mat 

x=0:L/numel:L; 
115        fori=l:length(x) 
h(i)=5; 

end 

clg 

hold  off 
120         %axis('off) 
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%plot([],[]) 
%hold  on 
plot(x,h,x,h,'x') 
125 
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%VoVo%%%%%%%%%%Vo%%%Vo%%%%%%%%%%%%%%%%%%%%%%%v0%%%%%%%%y0yo%%%%%%%% 
%%%%%%% 

%FORM  FE  BEAM  MODEL  WITH  ERRORS  AT  PRESCRIBED  LOCATIONS%%%%%%%%%%% 
%MASS  ERROR  LOCATION  SPECIFIED  VIA  POSOFMASSER.  VALUE  OF  ERROR  IN%% 
5        %PERCENT  GrVEN  BY  VALOFMASSERR,  ETC.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%o/o%%o/o%%%%o/o%%o/0o/0o/0o/0 

%%%%%%% 

clear,clc;clg 

casenanie='Spatially  Incomplete  A  set  Stiffness  error"; 
10        ifexist('setup.mat') 

load  setup  %load  existing  input  data  file 

else 
setup  %create  new  input  data  file 

end 
1 5        slruc_damping=0.0000000000 1 ; 

lumpdamp=luiupspring*0 

posofMasserr=[3  4]    ; 

valueofMasserr=[0.25  0.25]     ; 

posotStifferr=[3  4]; 
20        valueofSufferr=[0.25  0.25]  ; 

poso£Dainperr=[3  4]; 

valueofDamperr=[0.25  0.25]    ; 

for  i=l  :length(posofMasserr) 
lumpmass(posoiMasserr(i))=valueofMasserr(i); 
25        end 

for  i=l  :lengtli(posofStifferr) 
lwnpsprmg(posofStiffenti))=valueofSuiFerr(i); 

end 

for  i=l  :lengtli(poso£Dainperr) 
30  luinpdamp(posofDamperr(i))=valueofDauiperr(i); 

end 

elen=L/numel; 

ifbc=='pp' 
doftokill=2; 
35        elseif  bc='cc' 
doftokill=4; 

elseif  bc=='ff 
doftokill=0; 

elseif  bc=-cl' 
40  doftokill=2; 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o0/o0/o0/o0/o%%0/o%0/o°/o 

%BUILD  THE  ELEMENTAL  MASS  AND  STIFFNESS  MATRICES%%%%%%% 
45        %FOR  A  2DOF/NODE  "F.E.  ■  STRUCTURES.%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%ELEMENTSTffFTSIESSMATRiX%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
ke=[12  6*elen-12  6*elen; 
50  6*elen  4*(elenA2)  -6*elen  2*(elenA2); 

-12  -6*elen  12  -6*elen; 
6*elen  2*(elenA2)  -6*elen  4*(elenA2)]; 

ke=(eeii/(elenA3)).*ke; 
55        g=386; 

eleinke=eeii/(elenA3); 

%ELEMENTMASSMATRICE%%0/o0/o%%%%%0/o0/o0/o0/o%0/o%%%%%%%0/o0/o%% 
me=[156       22*elen         54         -13*elen; 
60  22*elen  4*(elenA2)      13*elen   -3*elenA2; 
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54     •   13*elen         156       -22*elen; 
-13*elen-3*(elenA2)     -22*elen  4*(elenA2)]; 

me  =  (((pho*area)*elen/g)/(420)).*me; 
65        elemnie=(((pho*area)*elen/g)/(420)); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o0/o0/o%0/o%%°/o% 

%%%%%%%%%%%% 

%ASSEMBLE  GOBAL  STIFFNESS  AND  MASS  MATRIX  FOR  A  2  DOF  F.E.  STRUCTURE%%%% 
70        %BASED  ON  THE  ELEMENTAL  MATRIXES  KE  AND  ME.  THE  STRUCTURE  CONSISTS  OF%%% 

%NUMEL  ELEMENTS  WITH  ELEMENT  CONNECTFVITY  GIVEN  BY  MATRIX  CONN%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o0/o0/o0/o0/o%%%% 

%%%%%%%%%%%% 
goblk=zeros(numdof); 
75        goblm=zeros(nuindof); 
goblc=zeros(numdof); 
for  i=l:nuniel 
v=conn(i,l); 
\v=coiin(i,2); 
80  goblk(dof_node*v-l  :dof_node*  v,dof_node*  v- 1  :dof_node*  v)=. . . 

goblk(dof_node*v-l:dof_node*v,dof_node*v-l:dof_node*v>i-... 
ke(  1  :dof_node,  1  :dof_node); 

goblk(dof_node*v-l  :dof_node*v,dof_node*w-l  :dof_node*w)=... 
goblk(dof_node*v-l  :dof_node*v,dof_node*w-l  :dof_node*w)f... 
85  ke(  1  :dof_node,dof_node+ 1 :2  *dof_node); 

goblk(dof_node*w-l  :dof_node*w,dof_node*v-l  :dof_node*v)=... 
goblk(dof_node*\v- 1  :dof_node*  w,dof_node*  v- 1  :dof_node*  v)+. . . 
ke(dof_iiode+l:2*dof_node,l:dof_node); 
90  goblk(dof_node*w-l  :dof_node*w,dof_node*w-l  :dof_node*vv)=... 

gob!k(dof_node*w-l  :dof_node*w,dof_node*w-l  :dof_node*w>+... 
ke(dof_node+l  :2*dof_node,dof_node+l  :2*dof_node); 

gobbn(dof_node*v-l  :dof_node*v,dof_node*v-l  :dof_node*v)=. . 
95  gobbn(dof_node*v-l:dof_node*v,dof_node*v-l:dof_node*v>4-.. 

me(  1  :dof_node,  1  :dof_node); 

gobbn(dof_node*v-l  :dof_node*v,dof_node*w-l  :dof_node*w)= 
goblm(dof_node*v-l  :dof_node*v,dof_node*w-l  :dof_node*w)+ 
me(  1  :dof_node,dof_node+ 1 :2*dof_node); 
1 00  gobbn(dof_node*\v-l  :dof_node*w,dof_node*v-l  :dof_node*v)= 

goblm(dof_node*w-l  :dof_node*w,dof_node*v-l  :dof_node*v)+ 
me(dof_node+l  :2*dof_node,  1  :dof_node); 
goblm(dof_node*w-l  :dof_node*w,dof_node*w-l  :dof_node*w)= 
goblm(dof_node*w-l  :dof_node*w,dof_node*w-l  :dof_node*w)+ 
1 05  me(dof_node+l  :2*dof_node,dof_node+l  :2*dof_node); 

end 

goblc=sqrt(- 1  )*struc_damping.  *goblk; 
goblkx=goblk; 
goblcx=goblc; 
110        gobbiix=gobbn; 

for  i=  1  :numel 
v=corui(i,l); 
115  \v=conn(i,2); 

goblkx(dof_node*v-l  :dof_node*v,dof_node*v-l  :dof_node*v)=... 
goblkx(dof_node*v-l  :dof_node*v,dof_node*v-l  :dof_node*v)+... 
lumpspring(i).  *ke(  1  :dof_node,  1  :dof_node); 

120  goblkxrdof_node*v-l:dof_iiode*v,dof_node*vv-l:dof_node*\v)=... 


140 


BEAMMDL.M 


goblkx(dof_node*v-l  :dof_node*v,dof_node*w-l  :dof_node*w)+-... 
lumpspring(i).*ke(l  :dof_node,dof_node+l  :2*dof_node); 

goblkx(dof_node*\v-l  :dof_node*w,dof_node*v-l  :dof_node*v)=... 
1 25  goblkx(dof_node*w-l  :dof_node*w,dof_node*v-l  :dof_node*v)+... 

lumpspring(i ).  *ke(dof_node+ 1 :2  *dof_node,  1  :dof_node); 

goblkx(dof_node*w-l  :dof_node*w,dof_node*w- 1  :dof_node*w)=... 
goblkx(dof_node*w-l:dof_node*w,dof_node*w-l:dof_node*w)+... 
1 30  lumpspring(i).  *ke(dof_node+ 1 :2  *dof_node,dof_node+ 1 :2  *dof_node); 

goblmx(dof_node*v-l  :dof_node*v,dof_node*v-l  :dof_node*v)=... 
goblmx(dof_node*v-l  :dof_node*v,dof_node*v-l  :dof_node*v)+... 
lumpmass(i).  *me(  1  :dof_node,  1  :dof_node); 
135 

goblmx(dof_node*  v- 1  :dof_node*  v,dof_node*  w- 1  :dof_node*w)=. . . 
goblnix(dof_node*v-l:dof_node*v,dof_node*w-l:dof_node*w>+... 
lumpniass(i).  *me(  1  :dof_node,dof_node+ 1 :2*dof_node); 

1 40  goblmx(dof_node*w- 1  :dof_node*w,dof_node*  v- 1  :dof_node*  v)=. . . 

gobbux(dof_node*w-l:dof_node*w,dof_uode*v-l:dof_node*v)+... 
lumpmass(i).*me(dof_node+l  :2*dof_node,l  :dof_node); 

goblmx(dof_node*w-l  :dof_node*w,dof_node*w-l  :dof_node*w)=... 
145  goblnix(dof_node*w-l:dof_node*w,dof_node*w-l:dof_node*w)+... 

luinpmass(i).  *rae(dof_iiode+ 1 :2  *dof_node,dof_uode+ 1 :2  *dof_node); 


goblcx(dof_node*v-l  :dof_node*v,dof_node*v-l  :dof_node*v)=... 
1 50  goblcx(dof_node*v- 1  :dof_node*  v,dof_node*v- 1  :dof_node*v)+ . . . 

luinpdamp(i)*sqrt(-l  )*struc_damping.  *ke(  1  :dof_node,  1  :dof_node); 

goblcx(dof_node*v-l:dof_node*v,dof_node*w-l:dof_node*w)=... 
goblcx(dof_node*  v- 1  :dof_node*  v,dof_node*  w- 1  :dof_node*  w)i-. . . 
1 55  lunipdainp(i)*sqrt(-l  )*struc_damping.*ke(  1  :dof_node,dof_node+l  :2*dof_node); 

goblcx(dof_node*w- 1  :dof_node*w,dof_node*v- 1  :dof_node*  v)=. . . 
goblcx(dof_node*w-l  :dof_node*w,dof_node*v-l  :dof_node*v)+... 
lumpdamp(i)*sqrt(-l  )*struc_damping.*ke(dof_node+l  :2*dof_node,l  :dof_node); 

goblcx(dof_node*  w- 1  :dof_node*  w,dof_node*w- 1  :dof_node*  w)=. . . 
goblcx(dof_node*w-l  :dof_node*w,dof_node*w-l  :dof_node*w)+... 

lumpdamp(i)*sqrt(-l  )*struc_damping.*ke(dof_node+l  :2*dof_node,dof_node+l  :2*dof_node); 
end 


nuindof=numdof-doflokill; 

170        [gk,gm,gc,k_aiial,m_anal,c_aiial]=fixbcs(goblk,goblm,goblc,bc); 

[gL\,gmx,gcx,k_exp,m_exp,c_exp]=fixbcs(goblkx,goblmx,goblcx,bc); 
save  beamdata 
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function  [kstat,mstat]=fstatic(k,m,oset,aset) 
aset_size=length(aset); 
kaa=k(aset,aset); 
kao=k(aset,oset); 
5        koo=k(oset,oset); 
koa=kao'; 
clear  k; 
k=[kaa,kao;koa,koo]; 

1 0        maa=m(aset,aset); 

mao=m(aset,oset); 

moo=m(oset,oset); 

moa=mao'; 

clear  m; 
1 5        m=(maa,mao;moa,moo]; 


20 


25 


t_static=-koo\koa; 
T_static=[eye(aset_size);  t_static]; 

kstat=T_static'*k*T_static; 
nistat=T_static'*m*T_static; 

end 
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% 

function  [kirs,nurs]=firs_tani(k,ni,oset,aset) 

% 

%  this  function  returns  the  IRS  reduced  stiffness 
5        %  and  mass  matrices,  given  the  unreduced  couterparts. 

%  Care  must  be  taken  that  the  aset  and  oset  vectors  correspond 

%  with  the  existing  arrangement  of  k  and  m. 

%  k  and  m  are  UNPARTTTIONED  matrices. 

% 
1 0        aset_size=length(aset); 

% 

kaa=k(aset,aset); 

kao=k(aset,oset); 

koo=k(oset,oset); 
1 5        koa=kao'; 

clear  k; 

k=[koo,koaJkaoJcaa]; 

% 

maa=m(aset,aset); 
20        mao=m(aset,oset); 

moo=m(oset,oset); 

moa=mao'; 

clear  m; 

m=[moo,moa;mao,maa]; 
25        % 

t_static=-koo\koa; 

Tstatic  =  [t_static;  eye(aset_size)]; 

% 

kstat=T_static'*k*T_static; 
30        mstat=T_static'*m*T_static; 

% 

tirs=t_static+inv(koo)*(moa+moo*t_static)*inv(nistat)*kstat; 

T_irs=[tirs;eye(aset_size)]; 

% 
35        kirs=T_irs'*k*T_irs; 

mirs=T_irs'*m*T_irs; 

% 

%  end  function  firs  tarn 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%%%%%%%% 

%%%%%%%%%%% 

%this  function  returns  a  vector  U  containing  modal%%%%%%%%%%%%%%%%%%%%% 
%frequencies  (rad/sec)A2  in  ascending  order  along%%%%%%%%%%%%%%%%%%%%%% 
5        %vviUi  associated  male  shapes%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0/0%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%% 

1 0        function  [u  1  ,lambda,index]=freqmode(k,m) 

[u,lambdat]=eig(m\k); 

[lainbda,index]=sort(diag(lainbdat)); 

ll=zeros(length(k),length(k)); 

for  i=l:length(k); 
15  ll(i,i)=lanibda(i); 

end 

lambda=diag(ll); 

forj  =  l:length(k) 
ul(:j)  =  u(:,index(j)); 
20        end 
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function  [r,c]  =  ndx3d(siz,ij,k) 

%NDX3D  Index  into  3-D  matrix  packed  in  a  2-D  matrix. 

%  ELEM  =  NDX3D([M  N  P],I,J,K)  returns  die  element  position  ELEM 

%  of  die  (ij,k)  elements  of  a  M-by-N-by-P  matrix  wliich  is 

5        %  stored  in  a  (M*N)-by-P  matrix.  For  example,  the  diree  M-by-N 

%  matrices  Al  ,A2 A3  are  packed  into  a  2-D  matrix  using 

%  A  =  [A1(:)A2(:)A3(:)]; 

%  If  lengdi(I)  is  m,  LENGTH(J)  is  n,  and  LENGTH(K)  is  p, 

%  Uien  ELEM  will  be  an  m-by-n-by-p  matrix. 
10        % 

%  [R,C]  =  NDX3D([M  N  P],I,J,K)  returns  die  row  and  column 

%  position  of  the  (ij,k)  element  as  stored  in  die  normal 

%  2-D  matrix  of  size  (M*N)-by-P. 

% 

15%  To  specify  all  die  elements  along  one  dimension  use ':'. 

%  For  instance,  NDX3D([3  5  4],2:3,':',3:4)  returns  die 

%  elements  for  die  2-by-5-by-2  matrix. 


20 


30 
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% 

%  See  also  ELEM3D,  MESHGRID,  SLICE. 

%  Clay  M.  Thompson  1 1  -3-92 

%  Copyright  (c)  1 992  by  The  MadiWorks,  Inc. 

%  SRevision:  1.7  $  SDate:  1993/09/03  14:36:52  $ 


25  if  isstr(i),  i  =  1  :siz(  1 );  end 
if  isstr(j),  j  =  l:siz(2);  end 
if  isstr(k),  k  =  l:siz(3);  end 


if  isempty(i)  |  isempty(j)  |  isempty(k),  r  =  [];  c  =  [];  return,  end 

if  any(  (i<0)  |  (i>siz(l )) ),  error(Tndex  I  out  of  range.');  end 
if  any(  (j<0)  |  (j>siz(2)) ),  errorChidex  J  out  of  range.');  end 
if  any(  (k<0)  |  (k>siz(3)) ),  errorCIndex  K  out  of  range.');  end 

35        ifnargout==2 

[jj,ii]  =  meshgrid(j,i); 
r  =  ii(:)  +  aJO-l)*siz(l); 
c  =  k(:); 

40        else 

[jj,ii4ck]  =  meshgrid(j,i,k); 

r  =  ii  +  (jj-l)*siz(l)  +  (kk-l)*prod(siz(l:2)); 


end 
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function  [K,M,C]=intsub(Z_I,  Z_R,  omega,  omega  1,  W,  setsize) 

waitbar_handle=waitbar(0,'Computing  Integrals'); 

j=sqrt(-l); 

5        integ=omega.A2.*abs(W); 
rint=real(integ); 
iint=imag(integ); 

aa=mytrapz(omega  1  ,rint)+-mytrapz(  omega  1  ,iint)*j ; 
integ=(  1  ./omega.  A2 ).  *abs(  W); 
1 0        rint=real(integ); 
iint=imag(integ); 

bb=mytrapz(omegal  ,rint>i-mytrapz(omegal  ,iint)*j; 
integ=abs(W); 
rint=real(integ); 
1 5        iint=imag(integ); 

cc=mytrapz(  omega  1  ,rint>i-mytrapz(  omega  1  ,iint)*j ; 
fori=l:set_size 
waitbar(i/set_size); 
for  k=i:set_size 
20  M(i^c)=(l/(aa*bb-ccA2))*(bb*mytrapz(omegal,... 

omega.  *Z_I(ndx3d([set_size  setsize  length(  omega )],... 
i,k,  1  :length(omega))).  *abs(W))-cc*mytrapz(  omega  1 ,. .. 
1. /omega. *Z_I(ndx3d([set_size  set_size  length(  omega)]... 
,i,k,l:length(omega))).*abs(W))) ; 
25  M(k,i)=M(i,k); 

C(i,k)=(  1  /cc)*mytrapz(omega  1  ,Z_R(. . . 
ndx3d([set_size  setsize  length(  omega )],... 
i,k,  1  :length(omega))).  *abs(W)); 
C(k,i)=C(i,k); 
30  K(i,k)=(l/(aa*bb-ccA2))*(cc*mytrapz(omegal,... 

omega.  *Z_I(ndx3d([set_size  set_size  length(omega)].. 
,i,k,l  :lengtl\(omega))).*abs(W))-aa*mytrapz(  omega  1 .... 
1. /omega. *Z_I(ndx3d([set_size  set_size  length( omega)]... 
,i,k,l  :length(omega))).*abs(W))) ; 
35  K(k,i)=K(i,k); 

end 
end 
close(waitbar_handle) 

40 
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function  z  =  trapz(x,y) 

%TRAPZ  Trapezoidal  numerical  integration. 

%  Z  =  TRAPZ(X,Y)  computes  tlie  integral  of  Y  with  respect  to  X  using 

%  trapezoidal  integration.  X  and  Y  must  be  vectors  of  the  same  length, 

5        %  or  X  must  be  a  column  vector  and  Y  a  matrix  with  as  many  rows  as  X. 

%  TRAPZ  computes  the  integral  of  each  column  of  Y  separately. 

%  The  resulting  Z  is  a  scalar  or  a  row  vector. 
% 

%  Z  =  TRAPZ( Y)  computes  the  trapezoidal  integral  of  Y  assuming  unit 

10        %  spacing  between  the  data  points.  To  compute  the  integral  for 

%  spacing  different  from  one,  multiply  Z  by  the  spacing  increment. 

% 

%  See  also  SUM,  CUMSUM. 

1 5        %       Clay  M.  Thompson,  1 0/1 6/90;  Cleve  Moler,  1/1 9/92. 
%       Copyright  (c)  1 984-94  by  The  Math  Works,  Inc. 

%    Make  sure  x  and  y  are  column  vectors,  or  y  is  a  matrix. 

20        %   Trapezoid  sum  computed  with  vector-matrix  multiply. 
[m,n]=size(x)  ; 
z=0  ; 

for  index=l:m 

yy=y(l+(mdex-l)*n:n+(mdex-l)*n); 
25        yy=yy(:); 

xx=x(index,:); 

xx  =  xx(:); 

z  =z+diff(xx)'  *(yy(l:n-l)  +  yy(2:n))/2; 
end 
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